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y  Acditional  Appendix  A> 

A  Two-Dimensional  Point-Mass  Program 
for  the  Passive  Receiver  Problem 

This  program  was  written  in  a  format  which  is  easily  adaptable  to 
an  arbitrary  two-dimensional  problem.  The  program  is  set  up  to  produce 
data  for  the  movie  of  Chapter  VI  .  Contour  maps  are  printed  to  provide 
information  regarding  the  density  shapes  and  data  is  written  on  magnetic 
tape. 

'  The  program  was  compiled  on  one  of  the  Kirtland  AFB  CDC  6600's 

using  the  RUN  compiler.  The  random  number  generator  used  is  simllai?’  fb- 
but  inferior  to  the  one  described  in  Chapter  IV.  !'  '  ■ 

Programmer;  K.D.  Senne 
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Adf’itional  Appendix  B. 


A  Gaussian-Sum  Program  for 
the  Passive  Receiver  Problem 

This  program  was  written  for  solving  the  "new  problem"  referred  to 
in  Appendix  B  of  Chapter  VI . 

The  program  was  compiled  on  one  of  the  Kirtland  AFB  CDC  6600 's  using 
the  RUN  compiler.  The  random  number  generator  is  identical  to  the  one 
used  in  the  program  of  Additional  Appendix  A. 

Programmer:  D.  L.  Alspach 

(Reproduced  by  permission) 
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Additional  Appendix  C. 

A  Gauss-Hermite  Program  for  Implementing 
the  Two-Dimensional  Phase  Demodulator 

This  program  was  written  to  test  the  Gauss-Hermite  method.  The 
numerical  results  of  this  program  are  given  in  Appendix  B  of  Chapter  VII. 

The  program  was  compiled  on  one  of  the  Kir t land  AFB  CDC  6600 ’s  using 
the  FTN  compiler.  The  code-optimization  algorithm  was  employed.  The 
randcm  number  generator  is  identical  to  the  one  described  in  Chapter  IV. 

Programmer:  C.  Hecht 
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Additional  Appendix  D. 

A  Point-Mass  Program  for  Implementing 
the  Interpolating  Version  of  the  Cyclic  Phase  Demodulator 

This  program  simulates  the  cyclic  interpolating  demodulator, 
described  in  Appendix  C  of  Chapter  VII. 

The  program  was  compiled  on  one  of  the  Kirtland  AFB  CDC  6600 *s 
using  the  FTN  compiler.  Tire  code-optimization  algorithm  was  employed. 
The  random  number  generator  is  identical  to  the  one  described  in 
Chapter  IV. 


Programmer:  C.  Hecht 
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Additional  ^pendlx  E. 

A  Fourier  Series  Implementation 
of  the  Cyclic  Phase  Demodulator 

The  Fourier  Series  experiment,  described  in  Appendix  D  of  Chapter  VII 
is  included  here. 

The  program  was  compiled  on  one  of  the  Kirtland  AFB  CDC  6600 's  using 
the  FTN  compiler.  The  code-optimization  algorithm  was  employed.  The 
random  number  generator  is  identical  to  the  one  described  in  Chapter  IV. 

Programmer:  C.  Hecht 
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Reprinted  from  Proc.  Second  S>mp.  on  Nonlinear  Estimation  Theory  and  Its 
Applications,  San  Diego,  Sept.  1971,  51-58. 


REALIZATION  OF  NON-LINEAR  FILTERS* 
R.  S,  Bucy** 


Abstract 

The  concept  of  numerically  exact  non-linear  filter  synthesis  is  defined  and 
shown  to  be  of  importance.  Confidence  intervals  for  error  performance  are 
described,  which  allow  one  to  make  meaningful  statistical  inferences  about 
filter  performance.  Further,  some  recent  synthesis  results  for  non-linear 
filters  are  discussed. 


1.  INTRODUCTION 

In  the  past  eleven  years,  a  rather  amazingly 
complete  theory  of  non-b’near  filtering  has 
emerged,  see  [ij  ,  [2],  [s] ,  and  [4].  Ii.  fact,  the 
representation  theorem,  see  [4],  has  provided  a 
solution  to  a  general  class  of  non-linear  filtering 
problems,  albeit  the  solution  depends  on  the 
evaluation  of  a  function  space  integral. 
Unfortunately,  the  applications  of  this  theory  have 
teen  almost  non-existent.  The  reason  for  this 
state  of  t’ffaii  s  is  paradoxical;  the  success  and 
wide-spread  use  of  the  Kalman-Bucy  linear  theory, 
see  [4]  and  [f].  The  situation  is  reminiscent  of 
♦he  popularity  of  the  frequency  domain  techniques 
in  the  mid-50  's,  wliich  delayed  the  state  soace 
viewpoint  for  some  years. 

Pn  this  paper,  I  would  like  to  review  some  of  the 
efforts  at  synthesis  of  optimal  filters  for  some 
specific  problems,  and  at  the  same  time  point  out 
some  slightly  paradoxical  things  I  have  learned 
over  the  past  three  years.  In  particular,  the 
meaning  of  Monte  Carlo  results  will  be  discussed 
in  some  detail,  as  it  is  now  commonplace  to 
design  systems  on  the  basis  of  Monte  Carlo  runs. 


too  often  only  25  to  100  Monte  Carlo  runs  are 
made.  Of  course,  as  we  shall  show,  such  a  small 
number  of  runs  produces  unbelievably  large 
uncertainties.  Finally,  we  will  discuss  in  detail 
a  parameter  estimation  problem  considered  by 
Lich*  in  [o]  and  the  phase  lock  study  in  [?]. 

2.  NUMERICALLY  EXACT 
REALIZATIONS 

The  basic  problem  of  machine  realization  of 
optimal  non-linear  filters  consists  in  the  iteration 
of  a  non-linear  convolution  equation  of  the  form 

Jn(xJ  =^Tn(x_,i) Jn-iU)  d^  (2. 1) 

with  *>n(x)  the  conditional  density  of  the  signal 
process  at  time  n  ,  given  the  observations  up  to 
n-1  and  Tj,  a  non-linear  function  of  the  last  piece 
of  data.  Of  course,  (2, 1)  is  nothing  but  a  discrete 
form  of  Bayes'  rule  which  can  be  written  as  (2, 1), 
when  the  signal  process  is  Markovian  and  the 
plant  noise  as  well  as  the  observation  noise  are 
independent  white  noise  sequences.  For  details 
on  how  can  be  rer  vsented  as  a  ‘S'^t  of 
(2  M+  weight!-  on  a  moving  grid,  see  [ll] 

and  [s].  The  important  point  is  that  the  point 


’’‘This  research  was  partially  supoorted  by  the  Air  Force  Office  of  ScienUf 
Research.  Office  of  Aerospace  Research,  United  States  Air  Force,  under 
Cf ant  AF-AFOSR- 1244-67. 

’•‘’’'University  of  Southern  California,  Los  Angeles. 

**■!'  d  is  the  state  dime. 
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mass  representation  is  numerically  exact,  that  is 
the  sup  of  the  difference  between  the  true  condi¬ 
tional  probability  of  a  compact  set  and  the 
approximating  probability  converges  to  zero  as 
the  number  of  point  masses  gets  large. 

To  illustrate  the  importance  of  numerical  exact¬ 
ness,  consider  the  following  problem- 


’‘o 


trial  and  the  statistic  is  computed  as 

o2  .  1 V  ~2 

^n  -  nx^i  ^i" 


with  X.  =  X.  -  X.  ,  then 
1  1  1 


, _  V  II  with  =  Ex*^  (3,1) 


'N  Vt'4-  f‘2 


"^n-1  ^  %-l  • 


with  u,  and  v.  independent  zero  mean  Gaussian 
white  noise  sequences  and  c  a  mutually  inde¬ 
pendent  zero  mean  Gaussian  white  noise  sequence. 
This  problem  was  considered  in  (ll)  and  (lOj. 

The  behavior  of  =  Ex^^l  z^^. . . . .  z^)  as 

computed  as  function  of  M,  the  number  of  weights 
in  tlie  approximation  to  i.s  quite  remarkable. 

As  long  as  none  of  the  observations  arc  large  in 
comparison  to  tlie  observation  noise  standard 
deviation,  an  M  of -It!  produces  x'j^y^(!VI)a»X|^yi^(2M) 
to  four  |)laces.  However,  if  z^  is  large,  then 
>{„.  (256)  =  .  What  happens  is  tl>at  in  the 

case  of  large,  the  filter  density  approximates 
a  delta  function  and  the  weight  approximation  to 
the  filter  density  must  contain  a  large  number  of 
weights  to  be  accurate.  The  phenomenon  is 
quite  interesting  as  it  shows  that  bandwidth  is 
irre'  evant  for  non-linear  filters.  Boc.ruse  we 
used  a  numerically  exact  realization  scheme,  we 
found  tliis  phenomenon,  which,  by  the  way  is  easy 

to  exnlain  after  it  was  observed  by  looking  at  the 
l/q 

error  of  the  estimate  z^^  . 

Without  a  numerically  exact  scheme,  there  is  no 
way  of  knowing  when  the  approximation  to  the 
optimal  filter  is  adecjuate.  We  note  in  passing 
that  the  sub-optimal  schemes,  such  as  |!2|  are 
not  numerically  exact,  a.s  is  re()laced  by  a 
Iruncated  se'  ies  approximation. 


is  asymptotically  normal  of  mean  zero  and 
variance  1  in  view  of  the  central  limit  theorem. 

p 

In  other  words,  S„  is  normal  mean  'nd 
^  n  d 

variance  ^  -1*2^) 


For  Monte  Carlo  error  analysis  then  is 
the  estimate  of  the  mean  square  error,  and  with 
probability  .  9972  _ 


-  ^n"^2-^ 


if  ^  3^2  as  it  i 
with  probability  ,  9972 

c  2 


as  it  is  when  x  is  Gaussian,  then 


H3Vh 


< 


3.  MONTE  CARLO  ERROR 
EVALl  A 'HON 

Suppose  that  a  rando;ii  \  iriabU  is  estimated  as 
X  and  that  \j  and  Xj  <ienot<>  the  i  independent 


In  Figure  I,  we  have  plotted  versus  N. 

Note  how  slowly  3  Vn  goes  to  zero.  The  net 
effect  is  that  in  order  to  acliieve  a  Zu  confidence 
interval  of  overall  length  .  2  of  ,  2,000 

Monte  Caro's  are  needed;  under  the  assumption 
u.  513^2^  .  Of  course,  if  p,j«3m,^,  the 

situation  improves,  but  if  (i,j  >  3ii^^  ,  the  situation 
is  worse. 

In  view  of  these  results,  the  Monte  Carlo  analysis 
of  systems  -which  is  often  made  with  IG  to  100  runs, 
seems  curious  at  best.  Of  course,  another 
advantage  of  synthesis  methods,  wliich  are 
numerically  exact  and  carry  the  distribution,  is 
that  M,j  can  be  at  least  appro.ximalety  deter¬ 
mined,  so  that  with  (3.  3)  the  confidence  of  an 
estimate  can  be  more  exactly  determined.  We 
used  2,000  Monte  Carlo  runs  to  evaluate  the 
performance  of  the  phase-lock  loop  in  (7)  wliile 
500  Monte  Car.o  runs  were  used  in  [fi] ,  as  w-e 
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proved  order  of  magnitude  mean  square  error  bet- 


and  if  such  a  graph  is  made  for  a  particular 
problem,  it  will  quickly  convince  one  that  no 
inference  should  be  made  on  ^*2  without  a  bare 
minimum  of  about  500  Monte  Carlo  runs. 


In  order  to  .ilustrate  that  the  problems  raised 
above  are  real,  we  consider  the  following  problem 
considered  first  by  Licht  in  [6] 


with  C  Gaussian  mean  one  variance  1  and  the 
dv  process  white  noise  of  spectral  density  R. 
Licht  used  the  Representation  Theorem  and 
obtained  the  experimental  results  shown  in 
Figures  II,  III  and  IV,  wliich  lie  disclosed  to  me 
in  [13],  after  I  had  seen  the  results  of  [s]  and 
suggested  he  simulate  the  sensor  orbit  filter. 
The  sensor  orbit  filter,  see  [4]  and  [o],  is  the 


the  optimal  versus  the  extended  Kalman- Bucy 
filter  as  the  confidence  intervals  are  too  large. 

In  fact,  none  of  the  experimental  results  are 

>^>{5 

precise  enough  to  draw  any  conclusions  . 
However,  the  true  sensor  orbit  performance  is 
significantly  liigher  than  the  experimental  extended 
Kalman-Bucy  filter  error,  on  the  basis  of 
l^gure  II  ,  only.  The  conclusion  that  we  may 
draw  then  is  the  following,  the  extended  K-B 
filter  performance  depends  on  the  coordination  of 
the  signal  process,  not  too  surprising.  However, 
even  more  significant  is  the  following: 

The  performance  of  the  extended  K,  B.  filter  may 
be  worse  in  sensor  orbit  coordinates.  (N,  B.  the 
wide-sense  sensor  orbit  filter  coincides  mth  the 
extended  K.  B.  filter  for  (3.  5)). 

Tills  result  surprised  me,  as  intuitively  I  expected 
the  extended  K.  B,  filter  to  perform  better  in 
sensor  orbit  coordinates,  I  feel  it  necessary  to 
remark  that  Lichl's  thesis  is  an  excellent  piece 
of  work;  its  omy  defect  consists  in  draivi  ig  con¬ 
clusions  from  too  small  a  number  of  Monte  Carlo 
runs.  This  defect  is  shared  by  many  papers,  for 
example  [l  2]  ,  where  50  Monte  Carlo's  are  used. 

3.  PHASE  LOCK  LOOP  DESIGN 


wide  sense  Kalman-Bucy  filter  for  the  system 


=  vyj  +  y2)^t 


where  y^  =  x,  y,  =  x^.  Note  that  (3.5)  is  just 
(3,4)  in  new  coordination.  The  true  error  per¬ 
formance  of  the  wide-sense  filter  for  (3.  5)  is 
easily  evaluated  by  solving  a  Riccati  equation  and 
is  shown  on  the  ligures.  Notice  first  if  we  assume 
that  'S  that  the  throe  plots  are  incon¬ 

clusive  cbout  the  relative  error  performance  of 


A  problem  which  has  an  interesting  liistory  as 
well  as  major  technological  importance  is  that  of 
phase  detection.  An  idealized  model  of  which  is 
the  followng- 

do  -  (Id  ;  Oq  ■  c 

dZj  -  cos  0  dt  +  dVj  (4.  1) 

dZg  =  sin  0  dt  4  dv, 

where  dd,  dv^,  dv,  are  mutually  independent 
wliite  noise  processes  of  spectral  densities  q,  2r, 
and  2r,  while  c  is  a  Gaussian  random  variable  of 
zero  mean  and  variance  A.  A  more  realistic 
model,  0  the  output  of  two  integrators  is  being 


■'  y«N(jJ,  P)  denotes  that  y  is  normally  distributed  mean  P  covariance  P. 
■  The  confidence  intervals  are  too  large. 
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studied  by  Hecht  in  [lO]  ,  howe  'er,  I  will  confine 
myself  to  a  discussion  of  the  model  (4. 1),  The 
interested  reader  should  consult  [?]  for  the 
detailed  analysis  of  the  optimal  filter  for  (4.1). 

It  can  easily  be  shown  that  the  steady-state 
extended  K-B  filter  for  the  model  (4. 1)  is  given 

by 


wi*>'  o(o)  =  0. 


do  = 


2  cos  0  +  dz^  sin  0)  ('',2) 


Now  0  given  by  (4.  2)  is  the  output  of  the  classical 
1  dimensional  phase-lock  loop  and  for  small 
^-0  the  error  is  R  =  V2rq  , 

Since  for  an  important  class  of  problems  only  the 
relative  phase  is  of  interest,  in  [v]  we  considered 
the  loss  function  L(0  -  0'  )  =  2(1  -  cos(0-0’*)) 
wliich  is  quadratic  for  small  o-  0'^  ,  but  cyclic 
otherwise,  in  other  words,  o'  was  chosen  so 
that 

E  L(0j^-<«.^*)lz.(j)  i  =  1,2;  j  =  l,..n-U 


<E  L(0„-0)(Zi(j)  i=l,2;  j  =  l,...n-l  ) 


(4.3) 


for  the  discrete  version  of  the  problem.  It  is 
easily  seen  that  is  given  by 


»  -1  ^n 

<'n  = 

'^n 


(4.4) 


with 


=  Esin0„|zi(j)  i=l,  2;  j  =  l, . . . ,  n-l) 
=  E  cos  0jjiz.(j)  i«^l,  2;  j=l, . . .,  n-l). 


Let  us  agree  to  call  the  cyclic  non-linear 
estimator.  It  is  easily  seen  that 


E  L(0„-0^')!z.(j)  i-1,2;  j=l  ....  n-l) 

==  2(1  -  ) 


(1.5) 


However,  since  the  lower  curve  represents  a 
lower  bound,  the  maximum  achievable  error  im¬ 
provement  would  be  2  dB.  The  optimal  filter 
achieves  1/3  to  2/3  of  the  error  performance 
difference  of  the  filter  which  senses  0  directly 
versus  the  classical  phase-lock  loop.  The 
maximum  achievable  error  improvement  for  the 
more  realistic  two-dimensional  phase-lock  loop 
(i.  e. ,  the  phase  is  doubly  integrated  white  noise) 
is  about  6  db,  so  that  one  might  hope  for  a  3  db 
error  betterment  for  the  optimal  filter  for  this 
problem,  and  this  would  certainly  affect  the  way 
on  how  phase-lock  loops  are  built  in  the  future. 
Unfortunately,  the  current  economic  doldrums 
have  adversely  affected  NASA's  ability  to  support 
such  a  study. 

4,  CONCLUSIONS 

In  the  past  three  years,  non-linoar  optimal  filters 
have  been  built,  using  digital  computers,  and  have 
been  evaluated,  using  Monte  Carlo  results.  The 
purpose  of  this  paper  has  been  to  indicate  the 
amount  of  Monte  Carlo  runs  necessary  to  make 
valid  statistical  inferences  as  to  filter  performance, 
as  well  as  reviewing  sor.  e  of  the  results  of  recent 
studies  of  non-linear  filters.  I  hope  this  paper 
will  have  the  effect  of  stimulating  meaningful 
work  in  the  area  of  non-linoar  filtering. 


In  Figure  V,  the  experimental  results  are  given 
together  with  the  theoretical  error  curve  derived 
by  Viterbi  in  [l4].  The  3  a  confidence  interval 
extends  .4  dB  above  and  below  the  non-linear 
error  performance  curve,  so  that  we  may  conclude 
that  the  optimal  filter  acliiever,  a  mean  square 
error  improvement  of  between  .  6  dB  and  1. 4  dB. 
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HYBRID  COMPUTER  SYNTHESIS 
of 

OPTIMAL  DISCRETE  NONLINEAR  FILTERS* 
R.  S.  Bucy,  M.  J.  Merritt  and  D.  S.  Miller** 


Abstract 

The  total  number  of  digital  operations  required  to  generate  an  estimate  using  a 
discrete  optimal  nonlinear  filter  is  extremely  large.  This  paper  delineates  the 
computer  task  and  introduces  an  algorithm  for  mechanizing  an  optimal  1-step 
predictor  using  a  hybrid  computer  to  reduce  the  computation  time.  The  advan¬ 
tages  and  limitations  of  the  hybrid  mechanization  are  described.  Hybrid  com¬ 
puter  predictions  (as  simulated  by  a  continuous  system  simulation  language)  are 
compared  to  those  generated  by  an  all  digital  mechanization.  Timing  estimates, 
accuracy  estimates  and  computer  equipment  requirements  are  given.  Exten¬ 
sions  of  the  algorithm  to  multidimensional  systems  and  observations  are  dis¬ 
cussed.  Use  of  a  parallel  digital  multiprocessor  to  perform  the  algorithm  is 
demonstrated.  The  use  of  interactive  computer  graphics  for  debugging  and  con¬ 
trol  of  the  algorithm  is  discussed. 


1.  INTRODUCTION 

It  has  been  known  for  a  long  time  that  the  optimal, 
discrete  time,  nonlinear  estimator  for  systems 
whose  plant  dynamics,  sensor  characteristics  and 
signal  statistics  are  known  is  just  the  algorithm 
obtained  using  Bayes'  Rule  to  sequentially  update 
the  conditional  probability  density  based  or.  the 
latest  data.  However,  until  recently  researchers 
have  avoided  this  numerically  exact  procedure, 
preferring  instead  to  develop  suboptimal  Taylor 
Series  estimators,  as  for  example  the  extended 
Kalman-Bucy  filter  and  the  minimum  least- 
squares  error  second  order  filter.  The  reason 
for  this  is  clear.  A  true  optimal  filter  is  a  com¬ 
putational  juggernaut.  In  the  past  few  years  the 
availability  of  3rd  generation  high  speed  digital 
processors  has  led  to  several  papers  by  Bucy  (3), 
Bucy  and  Senne  (4),  (5),  and  Bucy,  Se-ne  and 


Geesey  (6)  which  describe  digital  computer  me¬ 
chanizations  for  the  required  iteration  of  the  con¬ 
ditional  probability  density  function.  Both  one 
and  two  state  variable  nonlinear  problems  with 
Gaussian  statistics  have  appeared  in  the  litera¬ 
ture.  It  is  quite  clear  that  an  extremely  simple 
4  state  variable  filtering  problem  with  a  real 
time  data  rate  of  1  sample  per  second  pushes 
contemporary  high  speed  sequential  processors  to 
the  limits  of  their  capabilities. 

The  hard  core  of  the  problem  in  computing  the 
{n+l)st  estimate  is  the  need  to  re-evaluate  the  d- 
dimensional  conditional  probability  density  func¬ 
tion  each  piece  of  data.  This  requires 

the  evaluation  of  a  d-fold  convolution  integral: 

-R^ 

'l)di  (1) 


v:..  -r 
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If  the  statistics  are  Gaussian,  is  an  ex¬ 
ponential  with  vector-matrix  argument,  is 

the  previous  conditional  density  function  evolved 
to  time  t  .  If  the  J  's  are  M-ly  discretized  in 
all  directions  then  and  are  each  de¬ 

fined  on  a  grids  containing  points.  The  com¬ 
putation  requires  evaluations  of  the  integrand 
for  each  of  the  points  in  the  domain  of 
Thus  one  iteration  of  the  algorithm  requires 
(M'*)  '  =  exponentiations  each  of  which  entails 
the  computation  of  a  multitermed  quadratic  form. 
The  integration  and  statistical  computations  which 
follow  the  convolution  consume  considerably  less 
computer  time  proportionally.  Even  for  low  di¬ 
mensionality  filters:  2  state  variable  (2-D)  plant 
with  linear  additive  plant  d  observation  noise 
sequences,  Bucy  and  Senne,  (4),  introduced  ro¬ 
tated  coordinates  and  ellipsoidal  mode  tracking  to 
achieve  a  real  time  data  rate  of  one  prediction  per 
two  seconds.  Suppose  the  conditional  density  at 
time  tj^f  J^d)#  and  the  corresponding  element  of 
the  discrete  data  sequence  z^,  are  available. 

The  computation  may  be  speeded  up  by  generating 
the  exponential  portions  of  the  integrand  both 
simultaneously  and  continuously.  This  is  easily 
done  on  an  analog  computer.  However,  the  ana¬ 
log  computer  is  notably  poor  at  storage  of  time 
varying  hinctions  of  several  variables.  Analog 
computers  are  not  well  adapted  to  aritnmotic 
computations,  the  solution  of  algebraic  equations, 
etc.  Hybrid  computers  we.'e  developed  to  over¬ 
come  these  problems  by  combining  the  digital 
computer's  arithmetic,  storage  and  control  func¬ 
tions  with  the  analog  computer's  ability  to  solve 
large  numbers  of  simultaneous  differential  equa¬ 
tions  rapidly.  As  might  be  expected,  hybrid  com¬ 
putation  introduces  new  problems.  Some  of  the 
most  important  of  these  are: 


(1)  How  to  partition  equations  between  ana¬ 
log  and  digital  computers, 

(2)  Development  of  analog  representations  of 
filter  variables, 

(3)  Setting  of  initial  conditions, 

(4)  Effects  of  limited  bandwidth  of  analog 
computer  components, 

(5)  Effects  of  digital  to  analog  and  analog  to 
digital  conversion  times, 

(6)  Limitations  imposed  by  maximum  digital 
computer  I/O  data  rates, 

(7)  Delays  caused  by  analog  computer  mode 
switching  and  reset  times, 

(8)  The  effect  of  limited  analog  computer 
dynamic  range. 

These  problems  are  considered  in  later  sections  of 
the  paper,  and  estimates  of  performance  of  a  high 
speed  hybrid  computer  (containing  a  contemporary 
analog  computer  and  a  3rd  generation  high  speed 
digital  processor)  are  given.* 

The  sections  below  will  describe: 

(1)  The  computer  task  in  optimal  filter  syn¬ 
thesis, 

(2)  Review  the  digital  computer  algorithm 
proposed  by  Bucy  (Ref.  3)  for  a  one  state 
variable  (1-D)  problem  with  linear  plant 
and  cubed  sensor, 

(3)  Describe  the  hybrid  algorithm,  illus¬ 
trating  its  application  to  the  same  cubed 
sensor  problem.  The  basic  theory  and 
flow  charts  are  presented. 

(4)  Some  of  the  details  which  lead  to  a  dis¬ 
cussion  of  the  hybrid  algorithm's  speed 


*In  mid-1971  a  contemporary  analog  computer  might  be  an  Electronu  Associates 
E.AI  680,  Applied  Dynamics  AD-4  .,r  Conicor  Ci- 5000 .  The  high  speed  digital 
processor  might  be  a  Control  Data  CDC  6600  or  IBM  360/85. 
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and  accuracy  are  described. 

(5)  The  latest  results  obtained  for  the  cubed 
sensor  problem  (simulating  the  hybrid 
algorithm  with  a  digital  continuous  sys¬ 
tem  simulation  language)  are  presented. 

(6)  We  will  briefly  discuss  some  of  the 
interesting  extensions  of  the  above:  2-D 
and  multi-D  problems,  the  effects  cf 
nonlinear  plant  on  grid  spacing  and  grid 
floating;  and  the  use  of  a  parallel  digital 
multiprocessor  such  as  the  ILIAC  IV  to 
implement  the  hybrid  algorithm. 

(7)  Preliminary  timing  estimates  for  1-D 
and  2-D  filters  far  selected  hardware 
configurations  are  given  and  compared 
against  similar  estimates  for  the  all 
digital  algorithm, 

(8)  The  use  of  interactive  computer  graphics 
for  display  and  debugging  of  optimal 
filters  and  their  density  functions  is  de¬ 
scribed  , 

(9)  Finally,  we  summarize  the  advantages 
and  disadvantages  of  the  hybrid  algorithm 
as  it  compares  to  the  purely  digital  opti¬ 
mal  filter  realization  and  conclude  with 

a  description  of  work  in  progress. 


Subsequent  sections  describe  the  hybrid 
mechanization. 

In  the  all  digital  algorithm,  for  computational  pur¬ 
poses,  the  conditional  probability  density  function 
is  represented  as  a  set  of  point  masses,  i.e.,  a 
set  of  M  dirac  delta  functions  appropriately  posi¬ 
tioned  and  weighted,  see  Figure  (1),  In  order  to 
obtain  the  weighting  to  be  assigned  to  each  point  in 
the  new  grid  space,  a  convolution  must  be  per¬ 
formed  over  every  point  carrying  positive  mass  in 
the  old  grid  space.  This  convolution  is,  by  far, 
the  most  time  consuming  portion  of  the  calculation 
and,  in  essence,  is  "THE  COMPUTER  TASK"  in 
optimal  filter  synthesis.  If  the  running  variables 
of  the  conditional  density  function  J^are  finely  dis¬ 
cretized  over  a  large  domain,  then  the  computa¬ 
tion  can  become  burdensome.  If  the  grid  is  al¬ 
lowed  to  "float"  at  the  beginning  of  each  iteration, 
it  need  not  be  of  such  hip^  resolution  nor  broad 
span  as  would  otherwise  be  required.  This  leads 
to  a  considerable  reduction  in  computer  workload. 
In  (3),  the  grid  for  is  centered  at  the  mean  of 

J  ,  X  ,  and  adjusted  to  span  66  or  86  to  either 
n  n  n  n 

side,  whore  6  is  the  standard  d'jviation  of  J  , 
n  n 

2.1  THE  ONE  DIMENSIONAL  CUBED  SENSOR 
PROBLEM* 

The  model: 


2.  THE  COMPUTER  TASK  IN 

OPTIMAL  FILTER  SYNTHESIS 

Optimal  discrete  nonlinear  filters  are  mechanized 
as  an  iterative  process  wherein  the  same  basic 
procedure  is  recomputed  for  each  new  piece  of 
data.  During  each  iteration,  a  new  conditional 
density  function  and  a  new  conditional  mean  are 
evolved.  Many  substeps  and  numerous  calcula¬ 
tions  are  involved.  For  convenience  and  compar¬ 
ative  purposes,  only  the  one-step  predictor  is 
considered  here.  This  section  reviews  the  all 
digital  algorithm  of  Bucy  (3)  and  its  application  to 
a  one -dimensional  (1-D)  cubed  sensor  problem. 


ax  ,  + 
n-1 


^n-1 


X 


0 


(2.1) 


z 


V 


n  n  n 

The  signal  process  [x  }  ,  is  a  dis- 

Crete  time,  indexed  set  of  scalar  valued  random 
variables  with  a>  0. 


The  stochastic  process  {u  }  ,  is  a 

set  of  independent,  normal,  identically  distributed 
scalar  valued  random  variables,  with  variance,  q. 
Hence  the  conditional  density  of  the  (n-^l)st  state,  y, 
given  the  n^*'  state  is: 


*A  more  detailed  development  of  this  problem  is  given  by  Bucy  in  (3). 
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The  random  variable  c  is  a  scalar  independent  of 
the  process  and  having  a  density 


fc(§)  = 

^  V2nk 


The  observation  process  fs  }  ,  is 

a  scalar  valued  data  sequence.  The  stochastic 

process  {v  )  ,  _  is  a  set  of  indepen- 

n  n-  I , . . ,  ,m, . , . ,  ^ 

dent,  normal,  identically  distributed  scalar  ran¬ 
dom  variables,  with  variance,  r,  and  independent 
of  the  u  process  and  the  random  variable  c.  Hence 

the  conditional  density  of  the  n^  observation  given 
the  state  is: 

1  _L  ,  p3.2 

I  ,v  -  =  ;=e'2r  '  (2.4) 

Based  on  this  description  of  the  signal  and  obser¬ 
vation  processes,  a  sequential  application  of 
Bayes'  Rule  leads  to  the  following  update  rela¬ 
tionship  for  the  probability  density  function  repre¬ 
senting  at  time  t^^^j  conditioned  on  data 

z  ,..., Zn  which  we  abbreviate  as  J  ,,(y). 
n  u  nri  ' 

T  i..\  o  /  „  2q  '  2r  n 


,(y)  S  / 


J„(§)d? 


where  =  means  equality  to  a  normalization  con¬ 
stant  dependent  on  ....  Zq. 

There  are  essentially  ‘wo  problems  associated 
with  implementing  an  optimal  nonlinear  estimator: 

(1)  It  is  necessary  to  represent  the  non- 
parametric  conditional  density  in  some 
practical  storage  media  -  the  represent¬ 
ation  problem. 

(2)  Evaluation  of  the  integral,  equation  (2 . 5) , 
-  the  convolution  problem. 


a  §^(i)*  Considering  the  components  of  J  as 
njnnegative  masses,  §j^(i)  is  the  point  in  R^ 
th 

which  carries  the  i‘  mass  J  [§  (i)]  . 

The  continuous  probability  density  function 
is  thus  discretely  represently  by  M  delta  functions 
of  weight  §„(>)]  located  at  §jj(i)»  i=l, . . .  ,M. 

Now  define 

A  (2.0) 

S3  ^  S„0) 

and  (2.5)  can  be  rewritten  as  a  discrete  density 
function: 

1  ,  -  ,2  1  ,  =3,2 

,  ,  ,0^  ■2c<yi’®V  '2r‘‘‘n'^  ' 

J„+l(yi'=Ee  -  ■>  e  J  (?) 


with  i  =  1 ,  . . . ,  M 


where  ^  •’n+l'yi^  -  '• 

For  an  initial  condition  from  (2.3), 

-J-  5^ 

. A,' 

The  gridding  Ijj(j)  is  centered  at  mean  and 
spans  standard  deviations  lo  either  side 

where  N  is  determined  empirically.  Hence,  we 


^n^l'**  '  ’n  ■  ^'’n  ^  (M-l)/2  ' 


with  1  s  i  s  M 


and  §j(i)  =  IqI')  with  (from  (2.3))  XQ  =  Oand 


^0  --Vk. 


2,2  REPRESENTATION  With  the  representation  (2 . 6)-(2 . 9) ,  the  convolu- 

Bucy  (3)  describes  the  probability  density  function  (2.5)  is  solved  for  Jj|^j(y)*  After  normaliza- 

as  M  ordered  pairs  {J^Kjj(i)],  §„(»)),  1  --  1,...,  M  an'i  variance  arc  computed. 

with:  ?^(i)  a  map  from  the  set  {1 ,2 , , . . ,  m]  to  the  2,3  CONVOLUTION 

reals;  J  [?  (i)l  is  the  corresponding  value  of  J 

n  n  II  Given  the  reorcsentation  (2 .  fd- (2 . 4) .  t2  51 


2.3  CONVOLUTION 


Given  the  representation  (2.6)-(2.9),  (2.5)  is 
solved  for  >n  M  sequential  integrations  each 
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consisting  of  M  integrand  evaluations,  where 
■^n+l^'”^  prior  to  normalization. 

2.4  ESTIMATE  COMPUTATION 


The  hybrid  algorithm  derives  its  speed  advantage 
by  using  the  analog  computer  to  perform  most  of 
this  computation  in  a  paraT'el  and  continuous  mcin- 
ner. 


After  normalization  to  obtain  =°n^P“tation 

«  ->  2 
of  its  mean,  ^nd  variance  are 

straightforward : 

M 

^n4l  "  Vn+l<yi> 

1-1 

°n+l  =  ^yfjn+l<yi>^-  <^n+l>^ 

1=1 

An  overall  flow  chart  of  the  digital  algoiithm  for 
this  problem  is  shown  in  figure  (2) .  The  inte¬ 
grand  evaluation  and  computation  of  estimates 
are  shown  in  figures  (3)  and  (4).  An  expanded 
flowchart,  figure  (5),  shows  the  integrand  eval¬ 
uation  for  a  2-D  optimal  1-step  predictor  and  is 
included  as  a  further  illustration  that  convolution, 
and  in  particular,  integrand  evaluation  is  ihe.  prob¬ 
lem  in  discrete  optimal  filtering. 

3.  THE  HYBRID  ALGORITHM 

The  preceding  discussion  has  pinpointed  the  major 
obstacle  to  high  speed  real  time  filtering.  This 
obstacle  is  the  need  to  calculate  the  conditional 


3.1  THEORY 

This  section  describes  the  mathematical  basis  for 
the  hybrid  algorithm  and  its  architectures.  Par¬ 
ticular  attention  is  given  to  the  problems  associ¬ 
ated  with:  (1)  partitioning  the  computational  tasks, 

(2)  analog  and  digital  representation  of  variables, 

(3)  analog  computer  techniques  for  exponential 
teneration,  and  (4)  convolution  on  the  analog  com¬ 
puter.  For  convenience,  two  assumptions  are 
made: 

(1)  All  analog  and  hybrid  components  are 
perfect,  i.e.,  noise,  drift  static  and 
dynamic  errors  do  not  occur.  The  effects 
of  those  errors  are  treated  in  Section  4 
below. 

(2)  Analog  computer  components  are  repre¬ 
sented  by  their  block  structural  contin¬ 
uous  system  simulation  language  counter- 
pa ’•ts,  which  do  not  invert  and  which  do 
not  require  amplitude  or  time  scaling. 

3.1.1  Partitioning  of  the  Computational  Tasks 


density  functions  and  their  corresponding  argu- 

ment  values  times  for  every  piece  of  new  data. 

For  example,  in  the  cub^d  sensor  problem  the 

2 

integrand  must  be  computed  M  times  during  the 
M  evaluation  of  the  discrete  convolution 


M  __1. 
..(Yt)?  £  e  2q 

*  ^  i-1 


which  is  a  discrete  approximation  to  the  continuous 


Before  the  computational  tasks  can  be  assigned  to 
either  the  digital  or  the  analog  computer,  a  num¬ 
ber  of  important  properties  of  equation  3.1  must 
be  examined: 


(1)  The  integrand  of  the  discrete-  convolution, 
3.1a,  may  be  partitioned  into  three  pieces 

'2a 

(a)  ^ 

1  ,  .3,2 

(b)  Zp(|.)  =  e  j 


(0 

the  first  of  which  depends  on  yj  while  the 
others  do  not. 
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(2)  A  similar  partitioning  may  be  made  for 
the  continuous  convolution,  equation 
3.1b: 


(a)  Y^,(y,§)=e 


-t!:  (y-a§)2 


(b)  Z^(§)  =  e 


1  ,  -3.2 

-2;<v5  > 


(c)  J^(?) 

Again,  c~ly  the  first  component  depends 


(3)  If  r  and  q  are  of  the  same  order  of  mag¬ 
nitude,  thi;n  the  exponential,  will 

vary  sharply  when  |§  -z^|  is  small. 

This  wdll  not  be  true  of  Y^(|)  if  ai  1. 

(4)  Jjj(l)  is  the  result  of  an  arithmetic  ,om- 
putation  at  the  conclusion  of  the  preceding 
iteration,  or  an  initial  condition.  Further 

is  stored  in  tabular  form  o.n  a  grid 
which  will  not,  in  general,  be  suitable  for 
the  next  it'rativa  cycle.  The  values  of 
may  b."*  obtained  by  interpolation 
(or  occasionally  extrapolation)  on  the 
tabular  data. 

These  observations  suggest  the  following  parti¬ 
tioning  of  tasks: 

(1)  Use  the  digital  computer  to  compute  and 

store  the  terms  I  for 

j  =  1,2,,,.,  M. 

(2)  Use  M  similar  analog  computer  circuits 

to  generate  the  M  functions  Y^(yj.§) 
where  1  =  )  ,2, . , .  ,M  in  parallel,  where 
I  becomes  a  continuous  variable  corre¬ 
sponding  to  {§j)j,,,2 . M' 

(3)  Transmit  from  the  digital 

computer  to  the  analog  computer  main¬ 
taining  the  correspondence  between 
and  §  . 


(4)  On  the  analog  computer  multiply  the  M 
Yc(yi,§)  terms  by  a  piecewise  continu¬ 
ous  function  corresponding  to 

Zd(V 

(5)  On  the  analog  computer  perform  a  paral¬ 
lel  continuous  integration  of  the  product 
obtained  in  (4)  to  obtain  M  unnormalized 
points  of  the  new  conditional  density  func¬ 
tion,  J“^,(y). 

(6)  When  §  (on  the  analog  computer)  has 
covered  an  interval  corresponding  to  the 
domain  of  Zp(Cj)  Jj^(§j),  halt  the  analog 
computation  and  read  the  outputs  of  the 
M  integrators,  described  in  (5)  above, 
into  the  digital  computer, 

(7)  In  the  digital  computer,  iiorinalize  and 
store  the  new  discrete  conditional  den¬ 
sity  function  Jjj^j(y)  and  compute  its 

mean  i  ,  ,  and  variance  0 

n+ 1  n+ 1 

For  the  case  of  Gaussian  plant  noise,  the  analog 
computer's  tasks  may  be  separated  into  three 
pieces:  exponential  generation,  integrand  forma¬ 
tion,  and  integration.  The  relationships  between 
these  three  pieces  are  shown  in  Figure  6. 

3. 1.2  Analog  Computer  Representation  of 
Variables 

The  representation  of  ss  a  set  o*'  discrete 

pairs  was  described  in  Section  2,  above,  and  is 
illustrated  in  Figure  1,  The  function  Zj^iSj)  is 
stored  similarly.  Y^,(y,§)  i.s,  of  course,  an  in¬ 
finite  set  of  functions,  that  is,  there  is  a  function 
e'cry  fixed  value  y  '  yj.  M  of  these 
functions,  Yj^(yj,§)j  ,  j  ^  are  generated 

as  continuous  functions  of  time  on  the  analog  com¬ 
puter.  Four  typical  Y^(yj,|)  tarms  are  shown  .n 
Figure  7.  In  order  to  form  an  integrand  corre¬ 
sponding  to  the  discrete  terms  oi  (3.  la) 

=  1.2 . M 

1  =  1 , 2, , . . , M  or  to  tile  continuous  terms  of  (3.1b) 
Y^(y.  §)  /ij-(§)J^(5)  the  digitally  stored  terms 
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must  appear  in  the  analog  computer  as  fijnctions 
of  time.  Define 


Di  =  {CZj3(§:)  •  „ 

(3.2) 

to  be  a  set  tf  discrete  pairs  stored  in  tabular 
form  in  the  digital  computer.  Dj  may  he  converted 
into  a  piecewise  continuous  function  of  t  by  intro¬ 
ducing 

,  M 

Zpc(‘)  =  j>Jn<^j*  C“(t-ti)-u(t-t2n  (3.3) 

where  tj  s  t  <  t2 


of  its  impulse  response: 

-Ts 


1  -  e 


and 


‘2  =  Wl) 


where  the  digital  representation's  running  variable 
Sj  is  mapped  into  the  analog  computer's  indepen¬ 
dent  variable  time,  t,  by  the  mapping 


^n 

§— >*t 


which  consists  of  a  scaling  and  translation.  That 
is; 


If  the  sampling  interval  T  (seconds/ sample)  is 
small  compared  to  the  period  of  the  highest  fre¬ 
quency  in  the  signals  being  processed,  then  it  can 
be  shown* 


r'  t o\  ^  ..  “T6/2 


Thus,  the  use  of  a  ZOH  reconstruction  at  the  di¬ 
gital-analog  interface  results  in  an  almost  pure 
delay  of  ^  T.  For  a  more  detailed  discussion  of 
this,  see  Miura  and  Iwata  [9].  Since  Zjj(§.)  is 
known  beforehand,  it  is  easy  to  translate  it  for¬ 
ward  by  this  amount,  cancelling  the  delav  and 
eliminating  a  possible  source  of  errors.  This 
translation  is  accomplished  by  representing  the 
set  of  discrete  products  on  the  analog  computer 
by  the  time  function 


M 

lpctt]=  ^ 


(3.5) 


with 


t  = 

f  (?)  =  b  (1  +  c  ) 
n  '  n'  n' 

(3.4a) 

U[t]  = 

u(t-t 

j)  -  u(t 

where 

where 

‘l" 

‘<‘2 

with 

- 

¥- 

‘’n 

=  t/[2N6  ,4  2No  ,/(M-1)' 

n- 1  n-  1 

(3.4b) 

‘2  = 

n' 

c 

n 

1 

c 

1 

1 

c 

<D 

It 

(3.4c) 

'I  1  = 

■ 

For  convenience  below,  the  subscripts  c-  b  arc 
c  will  be  omitted,  and  it  will  be  understood  that 
b  and  c  are  functions  of  the  step  number,  n.  A 

I 

typical  Zp^  (t)  resulting  from  this  mapping  is 
seen  in  Figure  8.  The  output  process  of  the  di¬ 
gital  computer  is  characterizable  as  a  zero  order 
hold  (ZOH)  process.  The  transfer  function  of  the 
ZOH  may  be  derived  from  the  Laplace  transform 


as  shown  in  Figure  9  with  the  map 
Y 

n 

? - >1  . 

That  is , 

t  ^  l  .jl?)  =  b(?  4  c) 


(3.6) 


*A  factor  of  1/T  has  been  eliminated  for  convenience.  The  original  sampling 
process  contained  a  constant,  T,  which  cancels  it  out. 
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with 

b 

c 

where  b  and  c  are  in  general,  functions  of  n. 

This  tranformation  considers  the  probability  mass 

represented  in  the  digital  computer  by  a  point 

mass  located  3i  §j  to  be  represented  on  the  analog 

computer  by  the  area  under  a  step  function  located 

on  an  interval  whose  midpoint  is  at  t.  =  f  (§.). 

J  n  J 

For  future  use  we  add 

Y'^  ?(t)  =  ^  -  c  (3,7) 

with  b  and  c  as  in  (3,6), 


2Na  , 


[n5_ 


.  No  , 

4.  n-  1  1 

■  1  ■  (M- 1)  ■  *n-  1 


and 


Wj(t)  [yi  -  a§(t)]^  (3.10) 

Generalized  intecration.  Generalized  integration 
is  a  name  given  to  a  method  often  used  lo  generate 
explicit  composite  functions  on  the  analog  com¬ 
puter  (10),  (11),  This  technique  requires: 

(1)  The  derivative  of  the  argument  with  re¬ 
spect  to  t  be  available  during  the  genera- 
ol  the  function, 

(2)  The  initial  condition  for  the  function  be 
supplied  for  t  =  Iq,  where  Iq  is  the  time 
at  which  function  generation  starts. 

Let  s(t)  be  the  output  of  an  analog  integration 
ele.ment: 


These  definitions  differ  from  those  given  previ¬ 
ously  (Eq,  3,4),  All  '.ubsequent  reference.®  to 
Y^,  b  and  c  will  refer  to  Eqs.  3,5  and  3,6,  As  a 
final  point,  we  note  that  the  analog  computer's 
total  run  time  T,  which  determines  the  constants 
in  (3,6),  is  chosen  as  the  smallest  value  for  which 
accurate  hybrid  computation  can  be  performed  by 
a  given  computer  system  for  a  particular  filter, 
see  Section  4  below. 


3,1,3  Exponential  Generation 


Having  determined  that  the  exponentiation  required 
to  develop  the  terms  of  Eq.  (3, 1)  will  occur  on 
the  analog  computer  and  that  a  suitable  mapping 
of  the  variable  §  into  time  t,  is  given  by  (3.6), 
we  turn  to  a  discussion  of  the  actual  generation; 


YcCyj,l(t)1  ^e  '  (3,8) 

with  1  =  1 , 2 , , . . ,  M 

where  5(t)  is  given  by  (3.7),  This  is  depicted  in 
Section  1  of  Figure  6,  For  convenience  we  define 


e  =  e  ^ 


(3,^0 


or  4^  -  s  ®  0,  The  solution  of  this 
dw  dw 

auxiliary  equation  is  the  desired  continuous  func¬ 
tion.  Since  M  of  these  functions  must  be  gener¬ 
ated  simultaneously,  and  all  will  depend  on  their 
arguments,  yj,  it  is  not  possible  to  esta'  hsh  a 
correspondence  betv/uen  w  and  t.  In  ead,  gener¬ 
alized  integration  is  involved  which  depends  on  the 
following  application  of  the  chain  rule  for  differen¬ 
tiation  of  composite  funi  tions 

if 

J  J 

Jf  is  available  in  the  analog  computer,  then 
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Application  to  the  state  update  exponentials .  Ap¬ 
plying  the  method  to  3 , 9 

^  dv^ 

'’'l'*'  ■  dt  d§  dt 


from  (3.7),  = -g 


This  can  be  expanded  to 


(3.14) 


As  mentioned  above,  b  and  c  are  functions  of  n, 
as  are  the  y^.  To  obtain  the  initial  condition 
W  [§(t-))] 

e  expand  (3.9)  using  (3.7)  with  t  =  tjj 

-  1  o  --  iq  1  b  ,3  ,5. 


For  the  case  tg  =  0  this  reduces  to 

WjC§(0)]  Cyi  • 

/»  S  <»  * 


(3.16) 


Hardware  constraints  often  preclude  starting  the 
generation  of  the  e  ^  terms  at  t^  =  0.  This  is 
discussed  in  Section  3.2.2.  The  derivative  terms 
W^(t) ,  I  =  1,2,  ..,,M  shown  in  Section  I  of  Figure 
6  are  also  generated  by  the  analog  computer.  The 
constants  necessary  for  production  of  these  der¬ 
ivative  terms,  as  well  as  the  initial  conditions, 

WjC?(to)] 

e  ,  are  sent  from  the  digital  computer 

prior  to  starting  the  analog  solution.  Sec  the  dis¬ 
cussion  in  Section  3.2  for  additional  implementa¬ 
tion  details. 

3.1.4  Integrand  Formation 

Returning  to  Figure  6,  we  note  that  the  M  terms 

1-  1  Z  M  5®*^**°"  ^ 

each  multiplied  in  Section  U  by  the  same  term. 

The  output  of  the  digital  to  analog  converter  (DAC) 
in  Figure  6  is  set  to  the  proper  amplitude  every 
6t  seconds  under  digital  computer  timing  control. 


The  change  in  representation  from  the  set  of  M 
discrete  pairs  Dj  stored  in  the  digital  computer  to 
the  piecewise  continuous  function  Zp^[t]  during 
a  run  on  the  analog  computer  is  based  on  the  map 
and  constants  b  and  c  which  are  derived  from 
the  conditional  mean  and  variance  estimates  of  the 
previous  iterative  cycle. 

Subsequent  sections  discuss  the  manner  in  which 
this  digital  to  analog  data  transmission  is  per¬ 
formed  and  some  of  the  problems  and  limitations 
associated  with  it.  Assuming  Zpj,(t)  of  proper 
form,  the  analog  computer  performs  M  parallel 
continuous  multiplications  as  shown  in  Section  11 
of  Figure  6. 

3.1.5  Integration 

The  last  step  required  by  the  convolution  of  (3,1) 
is  the  M  parallel  integrations  of 

Wi(t) 

=12  M*  integrations 

proceed  as  indicated  in  Section  III  of  Figure  6, 
continuously,  in  parallel  and  simultaneous  with 
the  generation  of  the  e  terms  and  the  injec¬ 
tion  of  Zp^(t)  which  (see  Table  I)  is  identical  to 

Jn|n<«- 

If  the  observation,  ,  is  too  far  from  x^,  or  r  is 
too  small,  severe  scaling  problems  will  occur  in 
outputing  Zpj.(t)  to  the  analog  computer.  These 
rather  unlikely  events  are  treated  in  Section  3.2.2 
below.  We  note  also  in  Section  III  of  Figure  6,  that 
in  addition  to 


J  ft 


Zp^it)  dt 


the  analog  computer  is  able  to  generate 


W 

yi®  ^ 

•'  0 

2.U  ,  ,  r  2 

yi''n+i<yi>  7o 


Zpc(t)  dt 


"Zp^(t)  dt 


for  I  =  1,2, 
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at  the  expense  of  some  additional  equipment.  The 
judgment  on  whether  to  do  this  depends  consider¬ 
ably  on  the  computing  power  of  the  digital  com¬ 
puter  and  on  the  analog  equipment  available,  see 
Section  for  additional  discussion. 

3.2  REALIZATION  OF  THE  HYBRID  ALGORITHM 

There  are  two  steps  required  to  produce  an  oper¬ 
ating  hybrid  optimal  filter:  (1)  definition  of  t^e 
equations  to  be  solved  by  each  computer,  and  (2) 
development  of  software  to  link  the  two  computers 
together  realistically.  In  the  process  of  parti¬ 
tioning,  the  first  requirement  was  satisfied.  The 
reader  fai  iliar  with  hybrid  computati.on  will  rec¬ 
ognize  the  myriad  of  details  which  mast  be  care¬ 
fully  treated  in  linking  an  analog  and  digital  com¬ 
puter  together, 

A  flow  chart  of  the  hybrid  algorithm  (for  the  cubed 
sensor  problem)  is  given  in  Figure  10.  The  in¬ 
formation  paths  within  the  hybrid  computer  sys¬ 
tems  are  shown  in  Figure  11.  The  analog  com¬ 
puter  diagrams  for  the  convolution  integral  are 
presented  in  Figures  12  and  13. 


estimated  mean,  spanned  by  the  n'**  grid 

(see  Figure  1)*;  and  Jq  -  the  assumed  initial  con¬ 
ditional  density  function.  The  zerotl.  grid  base  is 
laid  out  using  Jq  as  a  reference.  Ti.'  observation 
number,  n,  is  set  to  -1.  The  processing  per¬ 
formed  here  is  identical  to  that  of  BOX  1  in  the  all 
digital  filter  algorithm  of  Figu'.  u  2.  The  digital 
computer  calculates  and  to  the  analog 

M  initial  values  of  yj  -  the  o  lin;  .-.t  which 
will  be  computed.  The  M  ini*..,  !  conditions 
W  (to) 

e  are  stored  in  ;he  an;.'-  ,  computer's  sample 

and  hold  elements,  as  are  b  :  .nitial  values  of  the 
scaling  constant,  bQ,  and  t) .  translation  constant 
Co.  Because  of  analog  cor  -  j  .tcr  limited  dynamic 
range  (see  Section  3,2.2),  Jine  initial  conditions 
will  be  too  small  to  be  adequately  set  on  the  analog 
computer  at  tii-ne  t  =  0,  '  n-:  digital  computer 
determines  the  time  in  ‘h  ■  olution,  Iq  ,  at  which 

each  of  the  M  exponential  generators  can  be  accu¬ 
rately  started.  When  t..is  time  is  reached,  the 
analog  gate  shown  in  Fr’ure  12  will  be  closed  and 
the  exponential  gener:iton  begun. 


3,2,1  Information  Flow  in  the  Hybrid  Filter 

The  hybrid  filter  is  best  described  by  progres¬ 
sing  through  the  flow  chart  of  Figure  10,  with 
occasional  references  to  the  flow  charts  of  the 
purely  digital  l-step  predictor.  Figures  2,  3  and 
4.  The  hybrid  flow  chart.  Figure  10,  assumes  a 
single  random  sequence  far  its  input.  In  general 
these  will  be  noisy  obsc  rvations  of  an  on-line 
stochastic  system.  The  functions  performed 
within  each  box  shown  in  Figure  10  are: 

BOX  1.  The  digital  computer  initializes  para¬ 
meters:  q  -  plant  noise  variance;  r  -  observa- 
noise  variance;  k  -  initial  state  estimate  variance; 
M  -  grid  discretization;  N  -  number  of  estimated 
standard  deviations,  0^-1’  (n-Dst 


BOX  2.  The  sample  ei^uence  number  n  is  incre¬ 
mented.  All  iterative  ;ycles  begin  here, 

BOX  3,  For  a  simulated  system  the  digital  com¬ 
puter  performs  a  model  and  density  update  iden¬ 
tical  to  that  reprsi  anted  by  BOX  3  of  Figure  2. 

In  the  ease  of  an  on-line  stochastic  input,  past  data 
corresponding  tc  the  incremented  count  is  assumed 
to  be  available.  However,  it  need  not  actually 
arrive  luitil  r<>  iJircd  by  the  computation  at  BOX  7. 

BOX  4.  The  analog  computer  is  placed  in  the  ini¬ 
tial  condition  (1C)  mode.  Initial  conditions  are 
set  into  the  analog  computer  integrators.  The 
analog  gates  in  Figure  12  are  closed  only  for  those 
integrators  at  which  tp  =  0.  Other  circuit  con¬ 
figurations  may  be  more  desirable.  The  important 


*As  indicated  in  the  discussion  at  BOX  5,  the  grid  base  is  not  changed  between 
iterations  unless  the  computation  requires  it.  Hence,  in  general,  N  is  the  num 
ber  of  estimated  standard  deviations,  about  the  (n-r)th  estimated  mean 

Xj^_^  spanned  by  the  nth  grid  where  r  is  the  number  of  times  since  the  grid  was 
last  changed. 
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point  is  that  limited  analog  computer  dynamic 
range  requires  that  come  integrators  start  "la.sjr" 
than  ethers  because  their  initial  value  does  not 
reach  an  accurately  settable  value  until  t^t^ 

where  t^  is  a  grid  base  and  problem  dependent 
number  (see  Section  3.2.2), 


BOX  5.  A  check  is  made  to  determine  whether 
the  grid  should  be  floated,  i.e.,  whether  the  mesh 
or  location  of  the  grid  base  on  which  the  (n+l)st 
probability  density  function  >s  to  be  con¬ 

structed  should  be  changed  from  that  on  which  the 
n*^  probability  density  function,  J^(§),was  placed. 
Although  much  time  can  be  saved  in  both  purely 
digital  and  hybrid  filtering  if  the  gj'id  base  does 
not  change,  the  percentage  savings  in  hybrid  com¬ 
putation  are  significantly  greater.  One  must,  of 
course,  ‘’e  very  careful  that  the  grid  mesh  is  not 
mislocated  or  too  na""  “ow,  causing  significant 
amounts  of  "probability  mass"  to  be  located  off 
the  end  of  the  grid.  Likewise,  a  grid  which  is  too 
course  will  not  be  able  to  properly  resolve  the 
probability  mass,  BOX  5  compares  with 

_  and  3  with  3  ,  where  r  is  the  number  of 

n-r  n  n-r 

observations  processed  since  the  last  grid  float, 
and  branches  to  BOX  6  if  either  of  the  differences 
is  too  great. 


BOX  6,  The  grid  is  floated  by  recomputing  the 
scale  factor,  b,  and  translation  factor,  c,  fol¬ 
lowed  by  computing  M  new  values  for  y,  and 
WjCtp)  I 

e  ,  This  information  is  DAC'd  to  tne  analog 


computer,  M  new  values  of  tg  ,  the  time  of  ap- 
th  ^ 

plication  of  the  1  initial  condition  e  to  the 
exponential  generator,  are  computed  and 


stored.  Interpolation  and  possibly  extrapolation 
are  required  for  because  values  of  |(j)  at 

which  its  M  constituent  delta  functions  ate  located 


should  correspond  with  the  new  values  for  yj.  For 
a  "smooth"  a  second  set  of  interpolations  is 

required  to  produce  (k- 1)  X  M  intermediate  values. 


The  k  XM  discrete  pairs  representing 


{JjjK(i)],  j  2  j^»  are  stored  for 

future  use  during  the  current  analog  solution  as 
well  as  for  all  future  solutions  on  the  same  grid 
base.  For  interpolated  values,  and 

define  a  set  of  discrete  pairs 


1,  2,  . 

{[Zj,(§,)V?,)],  I,},..  2, 


,,,  kXM 

(3.2a) 

kxM 


stored  in  the  digital  computer  corresponding  to 
Oj  (see  £q.  3.2).  If  k  is  chosen  odd,  the 
mapping  (3.6)  may  be  maintained  and  the  time 
function  appearing  on  the  analog  will  be 


I  S  1 


where  U,  is  defined  as  in  3.5,  but  with  running 
*  (1) 

subscript  I .  As  Zp^(t)  is  identical  in  form  to 
Zp^(t)  shown  in  Figure  9  (b\it  with  increased  dis¬ 
cretization),  the  superscript  will  be  suppressed  in 
this  sequel.  For  this  kind  of  intermittent  grid 
float,  the  span  and  center  of  ?  and  y  are  made 
identical  in  order  that  each  succeeding  J„^j(y)  is 
developed  on  a  grid  base  suitable  for  its  injection 
as  Jjj(?)  on  the  next  iteration. 

BOX  7.  For  an  on-line  system  there  may  be  a 
wait  for  the  n*^  data  sample,  z^.  Using  z^, 

Zp(?j)  is  computed.  The  product  Zj^(§j)  . 
is  formed.  Under  certain  conditions  (see  Section 
3,2.2)  those  values,  for  /.  >  1  (corresponding  to 
t  >  0),  may  be  computed  in  parallel  with  the  oper¬ 
ations  performed  in  BOXES  8-14  as  long  as  they 
are  ready  before  the  analog  computer  requires 
them. 


BOX  8.  The  run.iing  index,  i,  of  the  piecewise 
continuous  function  Zp^(t)  (see  Eq.  (3.5a)),  is  set 
to  1.  Zp(§j)Jj^(?j), or  what  is  equivalent ,J^jj^(§j) 
is  DAC'd  from  the  digital  to  the  analog  computer. 


69 


486 


BOX  9.  The  analog  computer  is  placed  in  the 
COMPUTE  mode.  This  is  time  tjj,  or,  as  is  the 
custom,  tp  =  0. 

BOX  10.  The  analog  computer  simultaneously 
and  continuously  generates 

e  =  e  ”  for  1  =  1  to  M>  where 

§(t)  -  ^  +  c,  see  Eq.  (3.7)  and  Figure  7.  In  addi¬ 
tion,  it  performs  a  continuous  integration  of  the 

W,(t)  .p, 

product  e  Zp^(t).  When  t  =  ]» 


i.e. ,  the  value  of  t  corresponding  to  the  time  of 
injection  of  the  next  digitally  stored  value  of 
Zpj,(t),  BOX  11  is  entered. 

BOX  1 1 .  A  test  is  made  to  determine  whether  the 

(k  X  value  of  J  |„(?,)  has  been  processed, 
n  jn  * 

BOX  12.  If  i<  k  X  M  the  index  t  is  incremented 
by  one. 

BOX  13.  The  next  digitally  stored  value  of 
■^nln^^^'  '^n|n^^/*  DAC'd  to  the  analog  com¬ 
puter. 


The  analog  computer  remains  in  COMPUTE  while 
I  is  incremented  from  1  to  k  X  M.  During  this 
time,  it  integrates 


operation  corresponds  to  the  one  taking  place  in 


BOX  4  of  Figure  2,  labeled  Integral  Evaluation. 
BOX  4  of  Figure  2  is  expanded  in  Figure  3.  It  is 
the  difference  in  nature  of  the  computation  being 


performed  by  the  hybrid  filter  in  BOXes  9-13  of 


Figure  12,  from  that  being  performed  by  the  dig¬ 
ital  filter  in  BOX  4  of  Figure  2,  that  leads  to  the 


major  differences  in  speed  and  accuracy  between 


the  two  optimal  filters. 


BOX  14.  If  the  answer  to  the  question  posed  in 
BOX  11  is  positive  indicating  that  k  X  M  time  sub¬ 
intervals  have  elapsed  (or  equivalently  that  t  =  T), 
the  analog  computer  is  put  into  HOLD  (i.e.,  inte¬ 
grator  inputs  are  disconnected). 

BOX  IS.  The  M  unnormalized  (n+l)st  probability 
density  function  points,  raad  from 

analog  to  digital  converters  (ADC'd), 

BOX  16.  The  digital  computer  performs  normal¬ 
ization  to  produce  )()')>  followed  by  computa¬ 
tion  of  X  .  ,  and  6  ,  This  calculation  is  identi- 
n+ 1  nr  1 

cal  to  that  performed  in  BOX  5  of  Figuie  2  and 
expanded  in  Figure  4. 

BOX  17.  A  check  is  made  to  see  if  the  last  piece 

of  data,  z  u  u  j  ^ 

n  has  been  processed.  If  n<nj^^j^, 

n  is  incremented  (BOX  2)  and  another  iteration  of 

the  algorithm  occurs.  If  n  *  •'mAX’ 

taken. 

3.2.2  Further  Details  of  the  Hybrid  Algorithm 

The  discussion  below  provides  additional  insight  in¬ 
to  the  functioning  of  the  hybrid  algorithm,  Analog 
and  hybrid  errors  are  treated  briefly. 

Grid  floating.  If  the  plant  is  unstable  (a>  1  in 
Eq.  2,1),  or  if  successive  estimates  do  not  con¬ 
verge  in  phase  space,  then  it  is  difficult  to  dis¬ 
cretize  the  conditional  density  function  “^^j(y)*  If 
the  range  on  y  is  made  large,  then  an  unreasonable 
number  of  discretization  intervals  results.  If  the 
discretization  is  made  coarser,  then  the  accuracy 
IS  impaired.  The  solution  to  the  impasse  is  to 
float  the  grid  (sec  References  3  and  4  ),  The 
floated  grid  permits  the  accuracy  of  a  finely  dis¬ 
cretized  grid  to  be  obtained  with  a  minimum  of 
computer  time,  by  placing  the  discretization  where 
it  will  do  the  most  good.  In  the  purely  digital 
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computation  if  the  grid  were  fixed  (i.e.,  not 
floated),  then  the  computations  described  above 
could  be  simplified  considerably.  The  functions 

<Vl  -  5j> 

1  ,  .3,2 

<='n* 

could  he  precomputed  and  stored  in  large  arrays. 
Some  interpolation  will  be  required  for  the  function 
which  contains  and  the  arrays  may  be  quite 
large.  Equivalently,  in  the  hybrid  algorithm,  the 
exponential  generators,  ?)>  could  be  estab¬ 

lished  for  the  entire  observation  sequence.  The 
sample  and  hold  elements  shown  at  the  output  of 
the  digital  to  analog  converters  (DAC's)  in  Figure 
13  would  not  be  necessary.  When  the  range  of  y 
in  J  ,  ,(y)  is  knovm  in  advance  and  is  not  too  large, 
acceptable  accuracy  has  been  obtained  using  both 
floated  and  nonfloate  i  grids.  As  the  range  of  y  in¬ 
creases,  the  difference  in  computation  time  for 
each  estimate  on  the  .  ited  and  nonfloated  grid  be¬ 
comes  larger,  at  equivalent  accuracy.  Hov/ever. 
those  problems  which  can  be  represented  by  a  cy¬ 
clical,  or  modulo-2n  density  function  not  requiring 
hi  ■'  local  resolution  such  as  the  phase  lock  loop 
....jr,  see  (12),  can  be  represented  on  a  fixed 
grid. 

Filter  errors  at  low  data  rates.  Analog  computer 
components  are  subject  to  a  variety  of  low  fre¬ 
quency  and  static  errors 

(a)  noise 

(b)  offset 

(c)  drift 

(d)  nonlinear  behavior 

Contemporary,  3rd  generation  analog  compatcrs(l9j 

are  available  with  low  frequency  and  static  at- 

curancies  on  the  order  of  0,01%.  As  the  number 

of  components  used  in  a  problem  increases,  these 

errors  arc  compounded.  Further,  the  range  of 

useful  computation  on  such  analog  computers  is 

_+  10  volts.  If  the  noise  level  is  approximately  1 

to  10  millivolts,  then  the  useable  dynamic  range 

4  3 

of  an  analog  computer  is  1  part  in  10  or  10  . 


This  limited  dynamic  range  should  be  contrasted 
with  1  part  in  10^^®  available  in  the  IBM  360 
series  machines,  and  the  even  greater  dynamic 
range  available  from  CDC  computers.  The  den¬ 
sity  function  of  a  normally  distributed  random  var¬ 
iable  drops  three  decades  in  3.2  standard  devia¬ 
tions  either  side  of  the  mean.  It  is  not  uncommon 
for  an  optimal  filter  to  span  12  or  16  standard  de¬ 
viations,  Digital  computers  easily  generate  densi¬ 
ty  functions  spanning  72  standard  deviations  36 
deviations  on  either  side  of  the  me.in).  Consider 
the  exponential  generator  in  Figure  6. 

To  generate  a  "bell”  on  the  analog  computer,  an 
'oitial  condition  is  needed  which  corresponds  to  a 
point  on  one  of  its  "skirts".  If  we  were  to  attempt 
to  set  all  initial  conditioti.s  at  the  beginning  of  an 
analog  computer  run  for  the  1-D  cubed  sensor 
problem  with  a  =  0.5,  fully  one  half  of  all  these 
initial  values  coi-respond  to  points  three  or  more 
standard  deviations  from  the  mean.  That  is,  half 
of  all  initial  conditions  would  have  to  be  set  to  an 
accuracy  below  the  analog  computer's  dynamic 
range.  To  avoid  large  errors  in  the  bells,  the 
start  of  integration  is  delayed  until  Y^(yj,?) 
reaches  a  value  v/hich  can  be  accurately  set  at  the 
integrator’s  output,  approximately  10  MV  m  a  10  V 
full  scale  .analog  computer.  This  value  is  stored 
in  the  tmple  .and  bold  device  until  the  time  corre¬ 
sponding  to  the  initial  condition  is  reached,  see 
Figure  12.  The  time  of  application  Iq  of  an  IC  is 
calculated  using  Eq.  3,  15  with 


W,[§(tQ)] 


=  ,001. 


th 

A  signal  Xp  (i)  from  the  I  '  digital  output  point  (DO) 
“I 

in  the  discrete  data  control  unit  (see  Figure  11)  re¬ 
leases  the  integrator.  If  all  of  the  decades  of  pre¬ 
cision  were  visible  in  a  graph,  then  the  exponential 
(unction  and  its  con esponding  analog  computer 
equivalent  might  appear  as  .follov/s; 


71 


488 


In  Section  5  we  report  results  of  a  hybrid  optimal 
filter  simulation  in  which  the  ratio  of  Vj^  to 
Vj^Ax  varied.  At  the  output  circuit  the  princi¬ 
ple  source  of  low  frequency  errors  is  located  at  the 
integrator  input  where  small  observation  noise 
variance  or  a  data  point  corresponding  to  an 

far  froit  x.  lead  to  a  needlclike  J_  i  (t).  Multi - 
Ti  n  |n 

pliers  are  built  such  that  they  produce  full  scale 
output  when  full  scale  input  is  applied  to  both  inputs. 

In  terms  of  threshold  voltages  (V_„),  if  the  full 
4 

scale  voltage  is  10  then  the  product  of  two 

inputs,  each  100  times  produces  only  a 

threshold  voltage  output.  Since  contemporary 
integrators  have  excellent  equivalent  input  drift 
and  noise  characteristics,  the  pi,mary  problem 
encountered  in  this  circuit  is  to  obtain  multiplier 
outputs  which  are  above  for  significant  inputs. 
Assuming  that  the  exponential  generators  are  al¬ 
ways  set  to  achieve  full  scale  voltage  at  a§  =  y^, 
only  the  scaling  or  method  of  multiplying  by 

J  I  (t)  may  be  varied,  to  avoid  loss  of  accuracy, 
n  In 

The  set  can  be  totally  precomputed  and  rescaled 
(at  BOX  7  of  Figure  10)  or  alternatively  multi¬ 
plying  DACs  instead  of  multipliers  can  be  used  for 
output  multiplication.  Both  slow  the  computation 
down  and  the  latter  increases  cost  as  well.  It  is 
also  possible  to  bypass  the  analog  computer  and 
perform  a  simple  analytical  computation  corre¬ 
sponding  to  a  stochastic  state  update  for  the  remote 
data  point  case.  See  (13)  for  a  discussion  of  the 


case  where  x  » 


1 

z  3. 


n  n 

Filter  Errors  Associated  with  High  Speed  Hybrid 
Computation.  Having  discussed  »o;.ie  of  the 


problems  which  arise  if  we  can  assume  the  analog 

computer  has  an  amplitude  accuracy  of  1  part  in 
3  4 

10  or  10  ,  it  is  important  to  determine  when  this 
assumption  can  be  made.  Above  a  certain  data 
rate,  finite  analog  computer  component  baudwidths 
and  nonzero  data  transmission  and  mode 
switching  times  will  so  deteriorate  analog  computer 
performance  as  to  invalidate  accuracy  analysis 
based  on  purely  static  errors.  Component  band¬ 
width  limits  the  fidelity  with  which  the  exponential 
generators  can  produce  Yj,(yj,?)  and  the  output 
multipliers  can  produce 

W 
e 

In  Section  4  we  present  analytical  .ip«ults  in  which 
the  exponential  integrator  is  replacci  Ly  a  second 
order  system  representing  Us  high  frequency  equi¬ 
valent  circuit.  The  output  multiplier  again  may  be 
approached  tnore  directly.  In  Section  4  we  present 
analytical  results  co’. sidering  two  limiting  cases 
for  the  form  of  j^(t) , 

Intercomputer  maximum  data  rates  and  hybrid  sys¬ 
tem  switching  times  are  a  second  source  of  errors 
in  high  frequency  hybrid  filtering.  The  hybrid 
optimal  filter  has  the  advantage  that  for  low  dimen¬ 
sional  problems,  i,e,,  d  s  3,  niay  be  com¬ 

puted  in  its  entirety  by  the  digital  computer  prior 
to  convolution  on  the  analog  computer.  That  is, 
the  two  computers  can  operate  essentially  sequen¬ 
tially,  This  can  be  seen  in  Figure  10  where  for 
n  a  0  the  digital  computer  performs  computations 
during  BOXes  16,  17,  2,  3,  5,  (6),  and  7  followed 
by  the  analog  performing  the  computation  in  BOX  10. 
The  linkage  system  and  the  digital  computer  do  not 
have  to  work  very  hard  while  the  analog  computer 
solves  simultaneous  differentia!  equations  at  high 
speed.  This  situation  does  not  prevail  in  closed 
loop  hybrid  computation.  .Maximum  channel  data 
transfer  rates  place  an  obvious  upper  bound  on 
filter  data  rate.  For  example,  to  produce 
on  the  analog  computer  in  1  MS  with  i  -  100 

J  I  (?j)  must  be  transmitted  at  a  100  KH  rate, 
n  |n  £ 

That  is,  every  10  ps,  a  word  located  in  computer 
memory  must  be  transmitted  to  a  register  at  an 
interface  (see  Figure  11)  and  conver  ed  to  an 
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equivalent  stable  analog  voltage  available  at  the 
DAC  output  shown  in  Figures  6  and  12.  Although 
higher  speed  equipment  is  available,  maximum 
channel  data  rates  appear  to  set  1  MS  as  a  lower 
bound  for  a  single  convolution  on  the  analog  if 
Zpj,(t)  is  of  the  form  shown  in  Figure  9. 

Even  though  the  analog  computer  performs  the 
convolution  in  an  essentially  open  loop  manner,  at 
high  repetition  rates,  phase  shifts  caused  by 
timing  imprecision  (i.e. , 'earliness  or  lateness) 
in  operation  of  integrators,  sample  and  hold  mode 
control  switches,  and  analog  gatis  are  an  additional 
cause  of  inaccuracy.*  Specifically  (see  Figures 
11,  12,  and  13),  these  errors  are  caused  by  the 
limits  on  timing  precision  that  the  X  (t)  signals 

can  close  the  analog  gates,  and  the  high  speed 
W^lt) 

DAC  can  supply  with  respect  to  each 

other  and  with  respect  to  the  time  origin  set  by 
the  COMPUTE  mode  control  switch  which  starts 
the  time  base  integ.'.itor .  We  muot  add  to  this 
timing  skew  and  the  effects  of  deteriorated  signal 
waveforms  of  the  pulses  transmitted  through  the 
linkage  system  to  initiate  these  operations.  It  can 
be  shown,  that  a  timing  imprecision  of  only  16 
nanoseconds  produces  a  .01%  F.S.  amplitude  error 
for  a  sine  wave  with  1  MS  period.  We  defer  de¬ 
tailed  examination  of  thin  point  until  Section  4  and 
mention  only  that  timing  imprecision  is  a  not  well 
advertised  hybrid  computer  characteristic  and 
that  a  300  nanosecond  switching  grey  band  is  a 
more  reasonable  value  to  expect  for  contemporary 
analog  computers  with  900  ns  -  1  Ms  total  switching 
times . 

Interpolation  of  As  a  final  detail  we  dis¬ 

cuss  the  advantages  and  difficulties  of  injected 
function  interpolation.  First  note  that  no  matter 
how  it  is  represented  in  the  digital  or  analog  com¬ 
puters  or  is  a  continuous  func¬ 

tion,  The  Y(~(yj,?)  terms  generated  on  the  analog 


are  also  continuous.  Therefore  interpolation  of 
Jj^(§j),  given  accurate  values  at  Ij  will  increase 
the  accuracy  of  computation.  Results  given  in 
Section  5  demonstrate  this.  However  interpolation 
lengthens  digital  computation  time  and  overlapping 
the  computation  leads  to  output  scaling  problems 
discussed  above.  In  addition,  as  mention,  the 
greater  the  number,  t  -  k  x  M,  the  more  data 
which  must  be  transferred  between  computers  per 
analog  run  and  the  larger  the  minimum  T  as  deter¬ 
mined  by  maximum  channel  data  rate  must  be. 
Another  approach  is  the  use  of  a  first  order  hold 
FOH.  This  device  (see  (16))  essentially  performs 
linear  interpolation  on  the  data  after  transmission 

It  does  not  require  a  “2  advance  and  requires 
less  transm  .tied  data  points  for  equivalent  ampli¬ 
tude  accuracy.  Only  a  simple  circuit  is  required. 
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4.  THE  EFFECTS  OF  HIGH  DATA  RATES 
ON  HYBRID  FILTER  ACCURACY 

As  mentioned,  there  is  some  data  rate,  abo''e  which, 
analog  computer  static  errors  are  no  longer  re¬ 
liable  guides  to  solution  accuracy.  In  this 
section  we  considei  for  the  scalar  filter  just 
what  the  tradeoffs  are.  Again  th:  cubed  sensor 
problem  serves  as  an  example. 

4.1  ANALOG  COMPUTER  COMPONENT  BANOWlOTH 

4.1.1  Exponential  Generator 

Under  certain  reasonabie  assumptions  an  operation 
amplifier  can  be  represented  by  the  following 


vFor  a  discussion  of  this  point  see  (14)  and  (15),  where  the  nature  and  effects  of 
H.  Y.  Brid's  law  also  called  the  "Magic  Equation"  arc  discussed. 
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transfer  function 

&U)  i  -  i!  . - ! - —  (•i.l) 

1  +■  -L- 

Aft 

where  A  °  open  loop  amplifier  Cm»i 

/Ui* 

feedback  network  short  circuit  transfer 
impedance 

il=  input  network  short  circuit  transfer 
impedrnce 

For  A  >>  I,  the  usual  approximation 

results.  At  high  frequencies  it  is  necessary  to 
replace  A  with  A(jjAa)  to  evaluate  operational 
amplifier  performance.  A(^xP)  is,  In  general,  a 
compi'  <  transfer  function  with  3  or  more  break 
points.  Simplified  analysis  valid  for  some  oper¬ 
ational  amplifiers  assumes  a  single  pole  transfer 
function 

Ao 

Substitution  of  <4.2  into  A.  I  and  setting 
/  It  \ 


(WHE-tl  'fcj  ^^.cRcl 

\  ^  ft  / 

for  an  integrator  yields  a  second  order  transfer 
function  containing  two  widely  separated  poles. 

The  low  frequency  pole  represents  the  physical 
fact  that  the  amplifier  gain  is  not  infinite  at  DC. 


f igurat ions 


rCn-O'+l 


bwiwrc  RPOCATC^ 


,  n  . 
i  JuJ 
II  !  i 


T'l  i\iM  I  nr  '  ■ 


t  (ortRPOLATEI) 


If  we  assume  the  fundamental  and  third  harmonic 
contai,-  all  terms  contributing  significant  mass  to 

manufacturer's  multi¬ 
plier  specifications,  a  A  hf  convolution  ti«e,  T, 
leads  to  errors  less  than  .02%.  Results  of  a 
simulation  of  this  portion  of  the  circuitry  will 
also  appear  in  (26) . 

4.2  HYBKID  SYbTEH  TIMIHG  IHIRECISION 
This  limitation  may  be  examined  By  considering  the 
mass  shift  in  ,  t']j)  ^Y  worst  case 

sketched  below.  1  ft-arl 

C*) 

pWiCt/  - L 


!^!|i  V 

0  A 


The  location  of  the  high  frequency  pole  determines 
the  effect  of  limited  bandwidth  on  hybrid  filter 
performance.  For  -  lo*  ,  >^,,0 

K)  "  «  ,  Tt  M  MS  cr~'  ^  ~  0 

It  can  be  shown  tnat  the  highest  frequency  of 
.0  oV/tC*) 

interest  In  \v  t  *•  will  not  produce  an 

error  in  excess  of  .01%.  Since  the  exponential 
generators  are  nonlinear,  tMs  analysis  can  only 
be  consiviereo  a  rough  guide.  See  (26)  for  results 
of  a  simulation  of  amplifier  limi'.sd  bandwidth  ir. 
an  exponential  generstO' . 

4.1.2  Output  Mu!:iplier 

_  Vj 

At  the  output  n'oltiplier  is  multiplied  by 


The  boxcar  function  corresponds  to  the  on  inter¬ 
polated  J^vj/n  v/orst  case  con,  luered  above; 

the  bell  is  ^  /j 

can  be  obtaince  rrum  norma! ited  tables.  For  the 
cooed  sensor  problem  with  1  »  4  AVs'l 

in  sketch  above,  'f  300  nanoseconds  (see  section 
3,2.2)we  again  find  rhat  the  errors  in  ) 

are  less  than  .( Z  1 

yS 

A  sir.ile.r  ^ort  of  dooiysis  In  wh'ch  an  Ideal 
J  Os  *  I  h  '■  coi  nred  against  the  JX, ,(<)).) 

produced  when  "X  %.)'s  delayed  bv  z;.j  jgs  aiso 


»  tv 

which  has  two  poss.bie  worst  ca:"  con-  yie'ds  erro's  in  -*”0  025 
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5.  SIMULATION  RESULTS 

A  simulation  of  the  hybrid  algorithm  for  optimal 
discrete  one  step  prediction  using  a  modified  ver¬ 
sion  of  MOE'SSL,  a  block  structured  continuous 
system  simulation  language  developed  at  the  Uni¬ 
versity  of  Southern  California  (17),  (18)  is  under¬ 
way.  The  modifications  to  MOBSSl<  provide  an 
ability  to  simulate  analog  computer  limited  dy¬ 
namic  range  and  make  double  precision  arithmetic 
available  for  highly  accurate  simulation.  Two 
types  of  results  are  av.iilable  at  present, 

5.1  V  ERIFICATION  OF  THE  ALGORITHM 

Table  2  compares  for  6  data  points  the  estimated 
means  and  variances  produced  by  a  simulated  hy¬ 
brid  filter  against  results  obtained  for  the  all  di¬ 
gital  filter  on  the  cubed  sensor  problem.  For 
both  filters:  k  =  q  =  r  =  1.0,  a  =  0.5  andN  =  8, 

For  the  hybrid  filter  M  r  ?9  and  k  x  M  =  319.  A 
dynamic  range  of  10  .i  v'as  used  in  this  simula¬ 
tion.  Digital  fil:er  .  discretizations  were 
M  =  29  and  M  =  j  Digital  filter  runs  of  the 
form  describo'i  bv  Bucy  (13)  to  obtain  Xjj(M)  'nd 
5^(M),  hav..  shown  that  for  this  problem,  M  =  513 
is  a  suftic'.en  ly  great  discretization  to  produce 
estimates  which  do  not  oiffer  in  the  4th  place  past 
ths  decimal  point.  Note  that  the  hybrid  filter  es¬ 
timates  are  closer  to  the  highly  discretized  digital 
f.lter  estimates  than  those  of  the  digital  filter  with 
M  =  29. 

5.2  EFFECTS  OF  ANALOG  COMPUTER  LIMITED 
DYNAMIC  RANGE 

Results  from  a  simulation  in  which  analog  com¬ 
puter  dynamic  range  was  varied  are  presented  in 
Table  3,  Mean  and  variance  estimates  obtained 
from  the  hybrid  filter  and  a  highly  discretized  all 

digital  filter  for  three  data  points  are  shown.  The 

2 

ratio  of  V., .  „  to  V_,,  is  varied  between  lO  and 
At  dynamic  ranges  of  10  or  greater  the 
‘imulated  hybrid  filter  estimates  are  in  close 
■  igreement  with  the  purely  digital  estimates.  For 
both  filters:  k  =  q=  r"l,a  =  0.5,N  =  8.  For 
the  all  digital  filters  M  =  513  and  M  -  29.  The 
hybrid  filter  grid  discretization  was  M  =  29  with 
no  interpolation,  i.e.,  k  X  M  =  29. 


6.  USE  OF  INTERACTIVE  GRAPHICS 
IN  OPTIMAL  FILTERING 

The  use  System  Simulation  Laboratory  contains 
a  Adage  AGT-10  graphic  computer  and  IBM- 360- 
44  medium  speed  scientific  computer  and  a 
Beckn  an  Ease  2132  Analog  Computer.  A  hybrid 
linkage  connects  the  computers.  Current  soft¬ 
ware  permits  the  graphic  computer  to  be  used 
standalone  or  in  conjunction  with  the  360  in  any  of 
several  modes , 

At  present  prC'grams  have  been  written  which  read 
a  sequence  of  stored  2-D  conditional  density  func¬ 
tion  from  disk  and  display  them  in  various  ways 
on  the  graphic  terminal  screen.  Most  of  these 
programs  are  written  in  the  AGNOS  graphic  pro¬ 
gramming  language  (20), 

These  display  programs  can  be  used  to  evaluate 
filter  performance  by  indicating  temporary  loss  of 
observation  data,  singularities  in  the  model  equa¬ 
tions,  grid  inadequacy,  poor  modeling,  or  very 
poor  prior  knowledge  of  state.  The  display  for 
one  of  these -programs  is  shown  in  Figure  14.  Any 
of  the  2-D  density  functions  on  the  periphery  may 
be  displayed  in  the  center,  A  sequence  of  mounds 
may  be  displayed  at  center,  using  a  "joystick"  to 
control  the  sequencing  rate.  The  density  functions 
may  be  blended  from  one  to  the  next.  The  graphic 
output  of  a  program  capable  of  displaying  any  pre¬ 
selected  density  function  is  given  in  Figure  15. 

The  graphic  computer  screen,  its  function  switches 
and  joystick  are  foreground.  The  IBM- 360  is  in 
background . 


75 


7.  SUMMARY  AND  CONCLUSIONS 
The  hybrid  algorithm  for  optimal  nonlinear  discrete 
filtering  has  been  presented.  Timing  estimates 
indicate  that  the  hybrid  algorithm  will  require  lO 
seconds  for  one  state  estimate  in  a  4-D  task.  This 
compares  to  3  seconds  for  49  CDC  6600's  operating 
in  parallel.  Thus,  one  hybrid  computer  with  250 
integrator;  and  multipliers  operates  at  speeds 
equivalent  to  49  CDC  6600's.  The  relative  disad¬ 
vantage  is  accuracy.  High  speed  hybrid  filter  es¬ 
timates  will  generally  be  within  1%. 

The  effects  of  analog  and  hybrid  errors  are  dis¬ 
cussed  and  a  more  thorough  treatment  is  in  pro¬ 
gress  [Reference  26]  . 

Application  areas  for  the  hybrid  optimal  filter  are: 

(1)  high  data  rate,  moderate  accuracy  tasks, 
i.e, ,  .  1%  to  1%  accuracy  is  acceptable, 

(2)  continued  Illumination  of  optimal  non¬ 
linear  filtering  with  a  view  towards  more 
general  theory, 

(3)  establishing  estimates  of  moderate  accu¬ 
racy  as  a  reference  for  suboptimal  filter 
estimates  in  high  dimensional  tasks. 
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11 

1.1 

1.2 

2. 

2.1 

3 

3.1 


TABLE  1  HYBRID  FILTER  NOTATION* 

The  Bayesian  Sequential  Estimator: 

=  Conditional  Density  of  random  variable  given  observations 

^0'  *1'  •••'  *b 

The  following  notational  simplification  is  used: 

If  a  =  b  +  conditional  subscript  b  ir  suppressed,  e .g. ,  ^  '^n+ 1  I  n^* ^ 

M  ' 

-'ntl<Vl>  =  2  Yj,(yj.|.)  Z^(g.)  J„(§.) 


Jn+i(y)  =r Yc(y-§) 

_ 1  <p 


1 

-jx  (y-a§) 

Y,.,(y,§)  =  e  ^  Unnormalioed  Conditional  Density  of  x  ,  ,  given  x 

^  11+  1  ^ 


-^(yi-a§)^ 

Y^.<y,,§)=e  2q  I 


Analog  Representation  of  1  for  fixed  y  =  y 


1  /  F  >2 

‘Iq 

YQ(yj,§j)  =  e  ^  •'  Discrete  Representation  of  1 

Zjj(§)  =  e  Unnorrnahzed  Conditional  Density  of  given  x^^ 


Zo(?j)  =  0 


Discrete  Representation  of  2 


J  (§)  Conditional  Density  of  Random  Variable  X  given  Observation  Zr., z z  ,  a  con¬ 
tinuous  function  ^  u  1  n-1 

Jn(§j)  Discrete  Representation  of  3 


4.1 


D  =  (iZ  (§  )  .  J  (5  )1,  §  } 

•'  J  J  ^  j  r  1  ,  2  .  .  .  ,  ,  ,V1 


f=l,2.  ...,kXM 


5.1 


Zpc(t)  =  ^  Zjj(§^)J^(5j)U(t)  Piecewise  Continuous  Function  Corresponding  to  4 
with  U(t)  =  u(t-tj)-u(t-t2)  '^here  tj^t  <  t2  tj  =  t2  = 


Y  defined  in  (3.6) 

II  \  • 


^PC  ^  Piecewise  Continuous  Function  Corresponding  to  4.1 


kXM 


6 

6.  1 


For  k  odd  define  U(t)  as  in  5  with  j  replaced  by  t 

^  Reference  8,  Eq.  2.5) 

■^nln*^''”  ^'D^^i*  "^n^^j'  D  screte  Representation  of  ( 


“'Unless  otherwise  specified  in  this  table  or  text,  the  iiota'ion  is  that  used  in  (3), 
(4)  and  (8). 


77 


494 


I-,-  -J-,;  <■;  - 


7 

W  (t)  . 

e  =  Y^(yj,§) 

Alternate  form  of  1.1 

TABLE  1  (CONriNue>) 

8 

W-U).  M 

jzl  J 

J„(t)  ^  E  u(t) 

j=i  ■' 

A 

Piecewise  Continuous 

Function  Corresponding  to  2.  1 

9 

Piecewise  Continuous 

Function  Corresponding  to  3. 1 

10 

U(t)  Piecewise  Continuous  Function  Corresponding  to  6,  I 

Note:  ZpcU) 

11 

•^n+l^yi’  '/q  ''^C^PC‘’‘ 

XH 

Unnormalized  !■*  Point  of  (n+l)st  Conditional  Density  Function 

Jn+l<y> 

TABLE  2 


COMPARISON  OF  ESTIMATES  FROM  TWO  ALL 
DIGITAL  OPTIMAL  FILTERS  WITH  THOSE  FROM 
SIMULATION  OF  AN  IDEAL  HYBRID  OPTIMAL 
FILTER 


ESTIMATED  MEANS  x. 


k+1 


ESTIMATED  VARIANCES  0 


k*  1 


DIGITAL 

HYBRID 

DIGITAL 

HYBRID 

M  =  29 

M  =  513 

M  =  29 

M  =  29 

M  =  513 

M  =  29 

0 

.855 

.825 

.825 

1.005 

1  .004 

1.004 

1 

-.026 

-.022 

-.022 

1.093 

1.092 

1.092 

2 

.334 

.333 

.333 

1.106 

1.097 

1.096 

3 

-.110 

-.100 

-.100 

1.099 

1.095 

1.095 

4 

.  145 

.  159 

.159 

1.089 

1.096 

1.095 

5 

-.012 

-.015 

-.015 

1.082 

1.085 

1.0S4 

TABLE  3 

EFFECT  OF  ANALOG  COMPUTER  LIMITED 

DYNAMIC  RANGE  ON  HYBRID  FILTER  ACCURACY 

ESTIMATED  MEANS  x 
DIGITAL 


k+I 
HYBRID 


ESTIMATED  VARIANCES  o 


•2 


k+I 


DIGITAL 


HYBRID 


M=513 

M=29 

10^:1 

10^:1 

lO'^il 

lO*'?! 

M=513 

M=29 

10^:1 

10^:1 

lO'*:! 

10>°:1 

0 

.825 

.855 

.845 

.854 

.844 

.854 

1.004 

1.001 

.869 

1.003 

1.0u9 

0.989 

1 

-.022 

-.026 

-.036 

-.039 

-.045 

-.020 

1.092 

1.0-, 

1.005 

1.090 

1.091 

1.091 

2 

.3  53 

.334 

.389 

.321 

.320 

.342 

1.097 

1.  106 

.979 

1.098 

1.096 

1.104 

Figure  5.  Integrand  Evaluation  in  a 
2-State  Variable  Optimal  I-Step  Predictor 


Mill’  inn 


Figure  7.  Exponentials  Generated  on  an  Analog  Computer  Corresponding  to  Four  Typical 
Terms,  each  of  which  Represents  the  Contribution  of  Model  State  Update  Equation  to  the  I Point 
of  the  (n+l)*‘  Conditional  Oensi ty  Function. 


O 


Figure  8.  Piecewise  Continuous  Function  Zfj^^Appeari.ng  at  Analog  Computer  as  a  Result  of 
Digital  to  Analog  Conversion  in  Accordance  with  the  Map  of  Equation  (3. A)  for  the  Case  H«.7. 
Dotted  Curve  Represents  a  Mapping  by  V;.,of  the  Funct  ion  2,  (f)J(nwhich  Coincides  wi  th  Z  (t)at 
the  H  Computed  Points.  <  i  « 
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Figure  11.  Data,  Timing,  and 
Control  Paths  in  a  Hybrid 
Computer  Realization  of  the 
Optimai  Discrete  Nonlinear 
Filter  Flowcharted  in  Figure  10. 
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£iail££_l|"  Analog  Computer  Components  in  I  Central 
Cwputational  Element  Group  of  a  Hybrid  Optimal  Discrete 
Filter  for  the  Cubed  Sensor  Problem. 
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DIGITAL  REALIZATION  OF  NON-LINEAR  FILTERS* 

Calvin  Hecht 
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Los  Angeles.  California 

Abstract 

A  method  is  presented  for  solving  the  discrete-time  non-linear  estimation 
problem  using  Gauss-Hermite  numerical  integration  and  Hermite  expansions 
of  the  conditional  density  function.  Numerical  results  are  given  to  demonstrate 
the  effectiveness  of  the  method. 


1.  INTRODUCTION 

Equations  have  been  developed  which  give  recur¬ 
sive  relations  for  the  density  function  of  the  state 
vector  conditioned  on  the  observations  (Bucy'2)  ). 
The  equations  are  general  and  allow  non-linear 
plant  and  observation  processes,  and  .lon-Gaussian 
plant  and  observation  noise  terms.  The  equations 
involve  integration  over  d-dimension  space,  and, 
except  for  d  =  1,  and  in  some  cases  2,  the 
equations  have  had  no  practical  solutions. 

Various  approaches  involving  approximation 
techniques  and  the  treatment  of  problems  less 
gei  ^'ral  than  the  complete  non-linear  problem 
have  been  proposed.  These  approaches  include 
schemes  for  approximating  the  non-Gaussian 
density  function,  methods  of  linearizing  the 
system  and  treatment  of  problems  with  a  non¬ 
changing  state  vector.  Although  the  methods 
have  had  different  levels  of  success,  a  goon 
generalized  method  has  not  been  produced. 

In  this  paper  a  2-dimensional  problem  is  solved 
using  a  technique  which  appears  capable  of  appli¬ 
cation  to  a  six-dimension  problem.  The  capa¬ 
bility  of  the  method  depends  on  the  extent  of  the 
non-linearities.  The  technique  can  be  used  for 
severe  non-linearities,  but  the  required  number 
of  points  for  numerical  integration  and  the  number 
of  terms  for  the  expansion  of  the  density 
function  become  large,  restricting  the  use  for 
higher  dimension  problems. 

2.  BAYES  LAW  COMPUTATION 

The  equations  for  the  conditional  density  of  the 
state  vector  conditioned  on  the  observation  pro¬ 
cess  using  Bayesian  estimation  were  derived  in 
Bucy''*'.  The  relationships  were  for  discrete 
time  processes.  The  equations  are  also  given 
in  a  slightly  different  form  (avoidan  ce  of  the 
pseudo-inverse)  in  Hecht'®' .  The  model,  notation, 
and  equations  as  given  in  the  second  reference 
are  given  in  Tables  2. 1  and  2.  2. 

The  relations  as  given  in  "  able  2.  2  are  generally 
very  difficult  to  solve,  e.ther  in  closed  form  or 


nuinerically.  Most  of  the  recent  effort  in  non- 
r.iear  estimation  has  been  toward  solving  the  re¬ 
currence  equations  for  specific  applications. 
Several  papers  at  the  Symposium  on  Non-Linear 
Estimation  Theory  of  1970  numeri¬ 

cal  techniques  or  other  schemes  to  obtain  practical 
solutions.  It  can  be  seen  from  equation  (2.4),  for 
example,  that  in  order  to  construct  a  density 
function  Jn/n(y)  defined  on  a  set  of  points  yj,  it 
would  be  required  to  integrate  over  the  entire 
d-dimension  space  R°,  If  the  prior  density 
function  were  defined  on  a  domain  which  con¬ 
sisted  of  k-points  for  each  component  of  x,  i 
would  be  equal  to  and  the  number  of  computa- 
t-ons  {o^  one  update  of  the  density  function  would 
be  (k®)  .  For  k  and  d  larger  than  approxi¬ 

mately  20  and  2,  respectively,  the  computation 
./ould  be  impossible  with  contemporary  computers. 

Bucy  and  Senne^^^have  devised  a  method  of  de¬ 
fining  the  function  and  performing  the  integration 
by  selecting  the  domain  in  a  special  way.  Pv  the 
efficient  selection  of  the  points  on  which  to  aefine 
the  density  function  they  were  able  to  solve  a  two- 
dimensional  problem  using  only  (15)^  points  and 
achieved  an  accuracy  equivalent  to  approximately 
(56)  points,  with  a  corresponding  decrease  in 
computation  time,  Bucy,  V  iler  and  Merritt 
have  proposed  reducing  computation  time  by  per¬ 
forming  the  integrations  in  parallel  on  an  analog 
computer  as  part  of  u  hybrid  computation  airange- 
ment.  Lo'°)  and  Alspach  and  Sorensom' '  re¬ 
present  the  den-sHy  by  a  sum  of  Gaussian  functions 
as  another  approach,  although  the  true  non-linear 
estimation  is  not  achieved  in  the  latter  paper 
because  only  a  linearized  system  is  considered. 
Kisner  (“)  and  others  represent  the  density 
function  with  Hermite  function  expansions  to 
solve  the  non-linear  problem,  but  also  fail  to 
demonstrate  the  Bayes  rule  recurrence  relation 
with  a  non-linear  sys.em  and  plant  driving  terms. 

In  the  following  development,  several  techniques 
are  employed  for  practical  solutions  to  the  equa¬ 
tions  of  Table  2.  2,  The  next  section  gives  the 
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analytic  details  of  a  two-dimension  example 
where  the  observation  process,  hn(x)  ,  is  non¬ 
linear  and  the  plant  is  linear.  Both  the  plant  and 
observation  noise  are  Gaussian,  By  making  a 
Hermite  expansion  of  the  non-linear  prior  density 
function,  the  a  posteriori  function  is  computed  in 
closed-form.  To  obtain  the  Hermite  coefficients 
for  the  nev/  density  function  requires  numerical 
integration,  but  the  Gauss-Hermite  numerical 
integration  formulae  are  utilized.  The  Gauss- 
Hermite  formula  has  been  demonstrated  in  Hecht 
to  be  eight  times  more  effective  than  rectangular 


integration  for  a  one-dimension  example.  That 
is,  equivalent  accuracy  was  achJeved  using  one- 
eighth  the  number  of  points  to  define  the  function 
to  be  integrated. 

The  technique  for  the  two-dimension  case  can  be 
extended  to  higher  dimension  examples.  In  the 
lugher  dimension  case,  computation  time  is 
reduced  by  performing  a  portion  of  the  integra¬ 
tions  analytically,  with  a  potential  of  i  ducing  the 
computation  time  for  a  d-dimension  problem  to 
the  equivalent  of  a  d/ 2-dimension  problem. 


Table  2. 1  -  Model  for  Bayes  Rule  Conditional  Density  Recursion  Formula 


’‘n  =  *n^’‘n-l^  %-l^’'n-l^  “n-1 


hn<V  ^  % 


Notation 


=  a  sequence  of  d-dimension  random  vectors 
=  time  index 

-  -  f. _ - od 


^n^^'n-l^  =  a  funcMon  from  R  to  R 

(7  ,(x  ,)  =  a  function  from  R*^  to  dvr  matrices 

n-i  n-i 

I  Uj^!  =  a  process  of  independent  r-dimension  random  vectors  with  density  Pu|^(w). 

c  =  a  d-dimension  random  vector  independent  of  the  u  process  and  having  density 

Pc(X). 

=  a  sequence  of  s-dimcnsion  random  vectors 
d  s 

h  (x  )  =  a  function  from  R  to  R 

n  n 

jv  j  =  a  process  of  independent  s-dimension  random  vectors  with  density  p^  (o), 

independent  of  c  and  the  u  process, 


Table  2,  2  -  Conditional  Uensity  Recursion  For.-oclae  for  Bayesian  Rstimation 
:currence  Equation  for  ('ondifio’ia'  J;*  nsitv 


Recurrence 

Relation 

- 

Jn-l/n-l<^>-*'^n/n<y’ 

''n/n<y’ 

Jn/n-l<y>  -*Jn/n<y> 

'’n/n^y> 

Jn/n-l<y>  ■*  ‘^n+l/n<y> 

‘^n+l/n< 

•fn/n<y)  ■*  Vl/n<y’ 

^n+lln^ 

Note:  The  index  n  has  been  omitted  in  some  places  for  notational  simplicity. 

Notation: 

Ja^b^y)  ■  conditional  density  of  the  random  vector  given  obser\  ation  vectors  .  .  Zj^ 

k(n)  =  a  normalizing  constant  to  make  .!(•)  a  density  function 

/.d./(  ♦)dx  =  d-  fold  integra.iun  with  respect  to  tlK- variables  x^,  components  of  the  vector  x. 
p^^(')  =  density  of  the  il-d;mensjon  random  vector  ^'n-I 
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3.  TWO-DIMENSIONAL  REALIZATION 


A  non-Unear  filter  was  synthesized  to  solve  a 
phase  angle  receiver/demodular  problem.  The 
model  for  the  system  consisted  of  position  and 
velocity.  The  receiver,  described  in  Viterbi  (13), 
and  Mallinckrodt,  Bucy,  Cheng dOl,  was  modeled 
as  observing  the  cosine  and  sine  of  the  position 
coordinate. 


The  equations  to  describe  the  signal  process,  or 
plant  equations,  were 


’'1 

*2 


(3.1) 


where  eu^  and  u2  were  uncorrelated  Gaussian 
white  noise  processes.  These  equations  were 
converted  t'  discrete  time  and  took  the  form 
Xi(n)  =  xi(n-l)  +  X2(n-l)A+ eui(n-l)i 
X2(n)  =  X2(n-l)  +  U2(n-l)A 

Hecht(®)  gives  a  summary  of  modeling  with 
discrete  time  systems  and  the  conversions  be¬ 
tween  di'^crete  and  continuous  time  models.  A  is 
the  time  between  sam|)les. 

Rewriting  equation  (3.  2), 


(3.  2) 


x(n  +  1)  =  Fx(n)  +  G  u(n) 


f - ‘1 

o 

1 _ 1 

G  = 

c  o' 
0  1_ 

fxi(n)' 

1  “  = 

rui(n 

Lx2(n), 

J 

Lu2(n) 

The  observation  process  was 
z(n)  =  H(xi(n),  X2(n))  +  v(n) 


(3.3) 


(3.4) 


The  density  functions,  using  the  notation  of 
Table  2.  2,  were 


(3.5) 


1 


(3.6) 


,,  E(c^  01*“*)  0 


E{uo  ) 


e^q  0 

0  q  I 


R  = 


rE[vj2]  0 

.  r  0 

L  0  E[v22]_ 

0  r  _ 

In  the  signal  process,  the  position  Is  a  deteimiai- 
stic  function  of  velocity,  which  is  accounted  for  by 
letting  e  0. 

The  filter  problem  was  solved  by  applying  p>qua- 


tions  (2.5)  and  (2.  7)  of  Table  2.  2: 

«0  CC 


•Vl/n<y)  =  mjJ  J 

-<1 

(3.7) 

•*n/n<y>  -m  Pv<^n  ’  Jn/n-l^^^ 

(3.8) 

In  the  expression  pgu(*  ^  equation  (3.  i 

if  ,V-'  If.'t 

c  -►  0  to  obtain  /  . 

( 

^  c^fiUyi-Xai^l-xijexp 

)'  2^' 
(.21 

where  _  ^y2~^2>  ){3.t;) 


6(  •  )  =  dirac  delta  function. 

Equation  (3.  9)  was  put  into  (3,  7), 

■^n+l/n‘y>  = 

with  c.  a  generalized  constant  which  includes 
k(n)  and  c^;  in  general,  we  let  cf  be  a  generalized- 
constant  which  may  include  contributions  from 
several  sources. 


7  te  integration  with  respec*  to  x,  was  performed 
by  taking  advantage  of  the  6-  function  to  obtain 

Vl/n<y>  = 

(3.10) , 

Vn{y2)  '  ^^2  expj--|p  [(zl-cosyj)2+(z2-&inyi)^J 

•^n/n-llyg) 

(3.11) 


The  pair  of  equations,  (3. 10)  and  (3.  ll)  give  the 
required  conditional  density  filter  update  relations. 
The  original  two-dimension  problem,  involving 
the  estimation  of  a  2-vector,  was  reduced  to  a 
single  integration  due  to  the  fact  that  only  one  of 
the  states  was  driven  by  a  random  !.'rocess. 

In  general,  when  the  signal  process  is  a  dynamic 
system  with  the  elements  of  the  state  vector  re¬ 
presenting  position  and  velocity,  the  covariance 
matrix  will  be  singular  since  the  only  driving 
terms  are  acceleration.  Thus  to  solve  a  problem 
in  three-dimension  physical  space  requires  a  six- 
dimension  state  vector,  but  the  problem  can  be 
reduced  to  a  triple  Integratiqn  in  the  same  manner 
as  above.  Bucy  end  Senne'^)  mention  the  computa¬ 
tional  simplification  resulting  from  the  singular 
covariance  matrix.  It  Is  the  above  simplification 
to  which  they  are  referring. 


The  integration  of  equation  (3. 10)  is  by  no  means 
an  easy  task.  Normally,  when  numerical  Integra- 
tlo  » is  applied,  a  posteriori  density  function 


j  ,  (y\]  is  computed  on  a  set  of  points  yi(i), 

2/ 

y2(i).  This  function  becomes  the  prior  density 
function  for  the  subsequent  iteration.  What  is 
needed  for  equation  (3. 10),  however,  is 
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VnL  ’'2  J  • 


Jn/nL  ^2  J  '  the  function  is  necSs-j} — , 

with  the  first  argument  different  thanthe  one  for 
which  the  function  was  computed. 

By  use  of  Hermiltpoljnomial  expansions  the  above 
difficulty  is  surmounted  and  an  explicit  formula 
was  derived  for  the  new  density  function  which 
involves  no  integrations. 

Because  of  space  limitations  in  this  paper  the 
detailed  steps  of  the  development  are  not  given. 

An  outline  of  the  steps  follows  with  equation(3. 1 2) 
the  final  result.  We  assume  a  set  of  quantities 
are  available  at  a  particular  time  which  complete¬ 
ly  characterizes  the  conditional  density  function, 
and  demonstrate  how  these  same  quantities  are 
computed  for  the  following  time  increment. 

1.  From  the  prior  stage  we  have  the  means, 
covariance  matrix,  characteristic  values  and 
vectors  of  the  covariance  matrix,  and  a  series 
expansion  of  the  conditional  density  in 
characteristic  vector  coordinates  in  terms  of 
Hermite  polynomials. 

2.  Series  form  of  the  density  is  put  into 
equation  (3. 10)  and  tlien  one  takes  the  Fourier 
transform  of  both  sides. 

3.  Rework  the  result  of  Step  2  to  take  the  in¬ 
verse  Fourier  transform  to  obtain  Jn+l/n  • 

4.  Multiply  Jn/n-l  by  the  observation  process, 
equation  (3. 1 1),  to  obtain 

5.  Compute  moments  for  means,  covariance 
matrix,  and  for  the  coefficients  of  the  new 
expansion. 


j  2A  ^2 

j  1 

\-l  zr 


^  (Zj-cos  y,) 
Mzj-stnyilJ  j 

<'^2  yr  Ti^ifi 


yj-  OiTjV,-  Tgyj) 


<Vi-T3 


(sign  Tj)  (sign  T2)* 


y2-(Ti  VI 


yi>)*”k 


(3.12) 


(^2"'^^^2yr’^i''i^'*^  f3^’^4yi  ■'''3' 


where  all  of  the  undefined  terms  are  functions  of 
the  quantities  listed  in  Step  i,  and  Cg  is  a 
normalizing  constant  to  make  a  density 

function. 

Jn/n  <  as  given  by  equation  (3. 12)  is  the  re¬ 
quired  equation  to  update  the  conditional  density 
function  for  each  sampling  time.  To  complete 
the  cycle  it  is  necessary  to  approximate  the 
density  function  with  a  new  Hermite  series,  as 
characterized  by  a  new  set  of  coefficients,  bj,j  . 
This  operation  is  shown  in  the  next  section. 

4.  HERMITE  APPROXIMATION  AND 
NUMERICAL  INTEGRATION 

Under  the  proper  conditions  a  function  of  two 
variables,  f(xiX2)  can  be  approximated  by  a 
series  expansion  of  Hermite  polynomials  of  the 

T .  ,  .  ^ 

f(xjX2)ft:  y(xi,X2)  =  c  e  ^  X 

r=o 


br(  IIj,(o'jXj)»  H  (02X3^  (4.1) 

v/ith  the  coefficients,  b^.^  ,  determined  by 

OlOp  f 

"ri  '  ~r - 1 - J  J  f(xi  X2)”r('"l  •''1  > 

2*^  r'  2^1"' 

Hf(  (>|X2)dXjdX2  (4.2) 


The  approximation  is  the  best  in  the  sense  that 
//e  e^  1  f(xjX2)  -  y(XjX2)|  dxjdx2 


mimmum 


By  letting  f(xjX2)  be  the  density  function  of  the 
previous  section  with 

9  .  n  .  9  2 


2  1  2  1 
"  rrr  .  "2  "  r~ 


1  ■  — 2"  2  ■  - r  >  <^i  >  1^9 

^  2oi‘^  •  2(72  ^  ^ 

the  characteristic  values  of  the  covariance 
matrix,  the  conditions  for  the  approximation  to 
be  valid  are  easily  satisfied  (Hildebrand  (7)). 


Let  f(xjXo)  be  a  density  function  we  wish  to 
approximtue,  and  let  mj  =  m2  =  0  in  equation 
(4.1).  Then 

'’l'^  r  V  .  .  ‘‘1"2 

*’00  ~T-  J  J  '“’‘i  '>>'2  = 


2  2  2  2 
y(XjX2)  =  —  e  e  (4.1 

2  2 

Let  <7.  and  02  be  the  variances  of  f(xjX2), 
and  let 

2  1 


155 


509 


Equation  (4.  5)  is 

y(x.X2)  =  2 


X  2 
2^ 


- 

2FT2 


(4.  6) 


By  taking  only  the  zero'th  term  of  the  expansion 
we  can  get  a  Gaussian  fit  to  the  true  density 
function. 

In  addition  to  the  above  assumptions,  except  let 
mj  and  m2  >  0,  let  ijj  and  1)2  be  the 

means  of  f{xjX2)  and  p,.  the  i  -  j'th  central 

moment.  ■' 

/•°° 

Pij  -  J  J  »)iMX2  -TJ2)^f(XjX2)  dXjdx2 

(4.7) 

The  product  of  the  Hormite  polynomials  when 
expanding  about  the  expected  values,  can  be 
written 

r  t 

Hj,(aj(xj-r,j))»ri|(a2(x2-5)2))  =  E  L  a.(x-) 

i=o  j=o 

aj(f)oi*  Of2'’(Xj  -  Tjj)* 

(X2  -  il2>^  <4,8) 

Equation  (4.  2)  is  then 

J  r.  J  X ! -00-00  1=0  j=o 

■  '^7'  •  ^  ^  a.(r)a.(f)aj*  <V2^  p. 

2*^  r!  2‘ £'7r  i  =  -  j  =  o  *  •>  *  2  ij 

^  (4. 9) 

The  coefficients  ajir)  aj(/)  in  (4,  9)  are 
functions  of  the  Hermitc  polyno...ial  coefficients, 
and  need  to  be  computed  only  one  time  in  advance 
and  stored. 

The  recurrence  formula  for  the  Hermite  poly¬ 
nomials  is  (Hildebrand  <71): 

K  (x)  =  2xH  (x)  -  iir'H.  .(x)  \ 

‘  *  (  (4,10) 

lyx)  =  1  H,  (x)  =  2x  > 

Designating  c^^^^  as  the  coefficient  of  x"  in 

H  (x),  a.(r)  a  ii)  in  equation  (4.  9)  is  deter- 
inuied  from  ^ 

(4.11) 


a.(r)  ajt£)  =  c^. 

IS  generated  from  a  re 
which  was  derived  from  (4. 10). 


was  generated  from  a  recurrence  formula 


Comparing  equation  (4.  i)  with  the  conditional 
de.iSity  function  in  eigenvector  coordinates  from 
the  previous  section,  we  see  the  approximating 
polynomial  expansion  can  be  accomplished  by 
letting: 

1  ^  .  I 


■'I 


“2^ 


’»2 


2  Ao 


and  using  (4.  9)  to  compute  the  coefficients  b^,^  . 

Gauss-Hermite  numerical  integration  was  used  to 
perform  the  integrations  for  the  moments  of 
equation  (4.  7).  Tlie  conditional  density  function 
was  given  by  equation  (3. 12), 

Gauss-Hermite  integration  is  performed  accord¬ 
ing  to  the  following  formula: 

2  2 


-X 


- 


f(x 


Xg)  dXj 


dx., 


n 

E 

i=l 


n 

E 


0,3  f(Xi., 


X2.) 


(4.12) 


with  xij  and  xoj  the  roots  of  H^fx) 
and  wj  ^nd  Wj  tWe  corresponding  weighting 
functions.  The  numbers  wj  and  xj  are 
tabulated  in  several  places  for  a  variety  of  values 
of  n;  Shao,  Chen,  Frank  <11/  for  example. 

The  a  posteriori  de.isity  function  of  equation 
(3. 12)  consists  of  a  product  of  an  exponential 
function  and  a  polynomial  function.  The  expo¬ 
nential  function  has  two  contributions,  the  non¬ 
linear  observation  part  with  the  trigonometrical 
functions  and  a  Gaussian  part. 

To  perform  the  numerical  integration  of  equation 

(3. 12)  effectively,  it  was  fcamd  desirable  to  con- 
,sider  two  situations  and  to  treat  them  differently. 

1)  When  the  contribution  of  the  observation  pro¬ 
cess  was  small  (r  large),  the  non-Gaussian 
exponential  was  combined  with  the  polynomial 
function  to  construct  f(xj,  X2)  of  equation  (4. 12). 
The  characteristic  values  and  vectors  of  the 
i-emaining  Gaussian  part  were  found  and  used  to 
formulate  the  exponential  part  of  the  integrand  of 

(4. 12) .  2)  When  the  contribution  of  the  observa¬ 
tion  process  was  significant  (r  small),  the  non- 
Gaussian  exponent  was  linearized.  The  linearized 
exponent  was  combined  with  the  other  Gaussian 
exponents  and  treated  in  the  same  way  as  in  1). 
The  difference  between  the  non-linear  and  linear 
exponent  was  combined  with  the  polynomial 
function  to  construct  f(xj,  Xg)  of  equation  (4. 12). 

5.  NUMER.CAL  RESULTS 

The  equations  of  the  previous  sections  were  pro¬ 
grammed  and  a  limited  number  of  runs  were 
mpde.  An  extended  Kalman- Bucy  filter  (extended 
K-.B  filter)  was  programmed  for  the  same  pro¬ 
blem  using  the  same  noise  sequence  for  compari¬ 
son  purposes. 

Ru(is  were  made  using  the  non-linear  filter  with 
thtf; coefficients,  b  ^  of  equation  (4. 1  2),  up  to 
b„g.  It  was  found  mat  equal  accuracy  could  be 
oDiiiined  by  eliminating  liigher  order  cross- 
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coefficents,  so  by  adjusting  mj  and  tn2  inequa¬ 
tion  (3. 12)  so  that  r+jt  <  8,  a  reduction  in 
running  time  without  loss  of  accuracy  was 
obtained. 

The  numerical  integration  was  performed  using 
10-point  and  20-point  integrations.  For  good 
signal-to-noise  ratios,  the  lower  number  of  points 
was  adequate,  whereas  with  poor  ratios  more 
accurate  integration  was  needed.  This  was  to  be 
expected  since  for  good  signal/noise  the  equations 
are  close  to  being  Gaussian  and  the  Gauss- 
Hermite  integration  formulae  are  almost  perfect. 
The  basic  independent  parameter  for  the  numeri¬ 
cal  trials  was  the  continuous  linear  model  error 
variance  of  xj,  designated  Pijio).  fiie  discrete 
linearized  variance  of  the  one-step  prediction  is 
designated  pii(A),  and  in  the  limit  pit(A)->pj  j(o). 
In  Hecht'®)  it  is  shown  that  for  A  sma)i  compared 
to  the  filter  time  constant,  F  , 


Pll(A)  ~  pjj(o) 


(1  +|-) 


All  trials  were  made  with  A/  F  =  0. 1. 

(131 

Viterbi'  '  shows  Pu(o)  is  the  inverse  of  the 
effective  signal/ noise  ratio,  so  Pii(o)  has  the 
double  meaning  for  the  continuous  linearized 


model; 


1)  = 
2)  Pjj(o)  = 


;[(xi  -^i)2] 


sTn 


s 


Trials  were  made  with  sample  path  of  100  points 
(10  filter  time  constants).  The  variance  of  the 
initial  state  was  set  equal  to  r  in  every  case.  The 
signal  term,  q,  was  set  toO.Ol  for  all  trials. 


For  very  small  pjjfo)  the  extended  K-B  filter 
and  the  non-linear  filter  performed  almost 
identically.  For  pii(o)  =,0001  the  sequence  of 
estimates  for  the  same  signal  and  noise 
sequences  agreed  to  within  6  decimal  places. 

As  the  parameter  PijCo)  was  increased,  the  non¬ 
linear  filler  was  also  tested  using  only  one  coef¬ 
ficient,  boo.  The  reason  for  this  is  that  one  would 
expect  better  performance  from  the  one-coeffi¬ 
cient  filter  than  the  extended  K-B  filter  due  to  the 
integration  of  the  non-Gaussian  term.  It  was  de¬ 
sired  to  determine  the  performance  of  the  nine- 
coefficient  filter  since  all  of  the  non-linear  in¬ 
formation  is  carried  by  the  coefficients  not  in¬ 
cluding  boo.  By  checking  both  the  one  and  nine- 
coefficient  filter,  the  ad^tional  improvement  of 
the  higher  terms  was  evaluated. 

For  pii(o)<  ,04  (S/N  =  13  db)  the  average  vari¬ 
ance  for  a  100  sample  run  was  almost  the  same 
for  the  extended  K-B,  the  one-,  and  the  nine-co- 
efficient  filters.  For  pii(o)  =  .lefS/N  =  8  db),  the 
one-coefficient  filter  shows  an  improvement  over 
the  extended  K-B,  and  the  nine-coefficient  filter 
shows  an  equal  improvement  over  the  one- 
coefficient  filter. 

For  Pii(o)  =  .  25  (S/N  =  6  db),  the  non-linear 
filter  appears  sub-^tantially  better  than  the 
extended  K-B.  These  results,  although  for  a 
very  small  sampling,  do  indicate  an  encouraging 
trend.  They  are  summarized  in  Table  5. 1. 


Table  5. 1  -  Average  Variance  for  Extended  Kalman- Bucy,  One-Coefficient,  and  Nine- 
Coefficient  Non-Linear  Filters. 

Sample  Path  =100  points 


Pll(o) 

Qd 

Rd 

A 

K-B 

1  -Coef, 

9-Coef. 

.01 

,00126 

0.05 

0.  126 

6.  6852  E-3 

7.0754  E-3 

7,0782  E-3 

.04 

.00200 

0.  2 

0.  200 

6.  1703  S-2 

6.  1681  E-2 

6.  1711  E-2 

,16 

.00317 

0.8 

0.317 

2.0403  E-1 

1. 9556  E-1 

1. 8349  E-1 

.25 

.00368 

1.  25 

0.368 

7.  2118  E-1 

2.  9438  E-1 

2,  7058  E-1 
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WITH  NONLINEAR  ESTIMATORS* 
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Abstract 

•» 

Elaborate  and  extensive  experiments  with  various  digital  Implementations  of 
discrete-time  nonlinear  estimators  reveals  some  Interesting  and  Important  charac¬ 
teristics  of  nonllneai  estimators.  The  experiments  Include  extensive  Monte  Carlo 
comparisons  of  alternative  techniques,  sample  path  comparisons  of  optimal  and 
approximate  estimators  and  associated  graphical  illustration  of  relevant  con¬ 
siderations  and  details. 


1.  INTRODUCTION 

In  the  proceedings  of  the  first  Symposium  on  Non¬ 
linear  Estimation  Theory,  Bucy  and  Senne  proposed 
a  digital  realization  for  nonlinear  estimators 
(henceforth  referred  to  as  the  Floating  Grid 
Method)'’^.  The  basis  for  the  Floating  Grid 
Method  is  a  numerical  approximation  to  the  general 
Bayes-optimal  solution  to  the  nonlinear  estimation 
problem^  Thus  the  Floating  Grid  Is  one  of  a 
variety  of  possible  applications  of  numerical 
techniques  to  calculate  the  Bayes  integral  recur¬ 
sion  formula**,^.  The  alternatives  to  the  numerical 
techniques,  which  Involve  only  density  representa¬ 
tion  and  Integral  approximations,  are  the  various 
analytical  techniques,  which  Involve  functional 
approximations  of  the  underlying  physical  problem 
equations  and/or  moment  series  expansions  of  the 


conditional  densities ^  Lomraon  analytical 
approximations  Include  linearized  (extended) 
Kalman-Bucy  filters,  quasi-moment  filters,  second 
order  techniques,  etc.  Ibe  obvious  advantages  of 
analytical  approximations  to  the  estimation  problem 
are  the  relative  simplicity  of  the  computations, 
compared  with  the  numerical  approximations.  On  the 
other  hand.  Improvements  to  analytical  techniques 
are  relatively  hard  to  obtain  since  generalizations 
come  at  the  expense  of  extensive  analysis  and/or 
engineering  experimentation,  whereas  the  accuracy 
of  a  numerical  approximation  can  be  directly  re¬ 
lated  to  computation  time,  with  little  or  no  modi¬ 
fication  to  the  computational  algorithm  Itself. 
Thus,  while  practical  applications  will  frequently 
demand  the  simplicity  of  analytical  approximations, 
the  answers  to  questions  relating  to  the  behavior 


*  The  views  expressed  herein  are  tho.“c  of  the  author  and  do  not  necessarily  reflect 
the  views  of  the  United  States  Air  Force  or  the  Department  of  Defense. 
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of  optimal  estimates  which  can  give  considerable 
insight  into  the  selection  of  practical  estimators 
may  best  be  obtained  by  studying  a  suitable  numer¬ 
ical  approximation. 

In  this  presentation  a  report  is  given  describing 
extensive  additional  experimentation  with  the 
Floating  Grid  Method  and  comparisons  with  ana¬ 
lytical  approximates  including  the  linearized 
filter  and  the  generalized  linearized  filter  of 
Alspach  and  Sorenson^.  All  techniques  have  been 
exposed  to  the  same  data  records  both  for  sample 
path  studies  and  elaborate  Hontc  Carlo  performance 
comparisons.  As  a  result  of  the  simulations,  many 
conclusions  have  been  made  concerning  the  behavior 
of  aptltaal  estimates,  the  "engineering  modifica¬ 
tions"  to  analytic  approximations  which  improve 
performance,  and  experimental  methods  relating  to 
evaluation  of  estimators.  Sample  path  studies  are 
exemplified  by  a  movie  illustrating  the  evolution 
of  the  conditional  probability  densities  for  a 
passive  tracking  problem.  Monte  Carlo  studies  arc 
given  to  compare  the  performance  of  the  Floating 
Grid  Method,  the  linearized  estimator,  the 
"taodlfied"  linearized  estimator,  and  the  gener¬ 
alized  linearized  (gaussian-sum)  estimator  on  a 
similar  tracking  problem.  The  results  illustrate 
conclusively  that  additional  efforts  on  sophisti¬ 
cated  numerical  techniques  will  result  In  con¬ 
siderable  additional  Insight  and  experience  with 
optimal  nonlinear  estimators. 

2.  PROBLKM  DESCRIPTIO.'J 

Consider  the  non-llncar  dynamical  system  described 
by 


X  "f  (x'  )+o  (x  )u 
— o  — n^-1  n-n-i-n-i^ 


x  “C 
o 


and  u^,  v^  are  zero-mean  white  Gaussian  sequences 
with  covariances  Q(n),  R(n) ,  respectively.  The 
initial  state  vector  £  is  distributed  on  accord¬ 
ing  to  the  density  J  (•).  If  the  conditional 

—X) 

density  of  x^,  given  meesurerents  z^,r j , . . . ,z^_ 


n-i 


is  denoted  by 

^(y)dy^rob(x^cy+dy  l^.z  j , 


|Z  ) . 
-n-1 


(2) 


Chen  an  application  of  the  discrete  representation 
theorem  gives  the  sequential  formula  formula  for 
obtaining  (y)  from  ^(y)  as  follows: 


(3) 


z  "h  (x  )+v 
—a  — n  -n  — n 


(1) 


with 


f  (•)  :  rW 

-n 

o  (•)  :  R'*-r'*xR’’ 


h  (.)  :  R'^-r’ 
— n 


whe^rp^^  is  a  normalizing  factor,  and  T^  depends  on 

(c)  and  v-f  (;)  in  accordance  with  the  assumed 
jm  -n  —  —  — n  — 

/ Causslan  density  functions^.  The  approximation 
problem  is  Involved  with  representing  the  condi¬ 
tional  densities  (2)  numerically  and  calculating 
Che  d-  dimensional  Integral  (1).  In  many  cases  of 
interest  the  conditional  moan  is  the  de,slrcd  esti¬ 
mate  of  X  . 

—a 

3.  PROPOSED  IMPLEIS’.NTATIOSS 

Some  early  attempts  to  Implement  approximate 
optimal  nonlinear  estimators  resulted  In  the  appli¬ 
cation  of  power  scries  expansions  of  the  system 
nonllncarl clc.s  and/or  moment  expansions  of  tie  '•cn- 
dltlonal  density  lunctios  ’  • •  .  Such  techr.icwes 
may  be  referred  to  by  the  term  "analvtital  approxi¬ 
mations",  since  the  underl>lrg  problem  Is  reorc- 
sented  by  a  simpler  analytse  iiscrlpticn.  ,'centlv 
some  attention  has  been  tumej  coward  tnc  ,  roolcm 
of  representing  the  Bayes  Imttgial  solution  to 
dlscrctc-tlcc  estimation  problems  by  a  nun.erlcal 
approximation  on  a  digital  computer.  Numerical 
techniques  consist  basically  of  two  Interrelated 
sub  problems  -  density  repre'^e-tatl  m  a-  B.'ive 
Integral  comp Jt.itlm.  For  1  r.s!ty  representation 
various  Investigators*''^  have  proposed  orfiugonel 
series  rcpresent.itlons ,  floating  point  grid  te-ro- 
scr.tatlcas' ’  * ,  and  suru)  c’’  go_ssl.in.s  .  Ihe  ic-'-cs 
Integral  cal^j’.allcn  Is  a-.uai'.;  tied  .ery  cle.e.v 
to  the  type  of  density  representation.  K'r  exanple, 
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orthogonal  series  representations  usually  lead  to 
soBe  fom  of  fourler  transfom  calculations  for 
the  Bayes  Integral,  either  explicitly**,  or 
Implicitly®.  Alternatively,  the  fleeting  point 
grid  representations  typically  lead  to  simple 
quadrature  formulas  for  the  Bayes  Integral.  To 
date,  however,  the  gausslan-sum  representation  has 
required  the  use  of  analytical  approximations  such 
as  local  linearization.  In  order  to  effectively 
implement  the  Bayes  law  computation,  which  really 
classifies  the  technique  as  a  generalized  linear¬ 
ized  (analytical)  technique. 

4.  AN  EXAMPLE  PROBLEM 


In  order  to  compare  the  performance  of  the  optimal 
estimator  with  that  of  currently  popular  approxi¬ 
mate  estimators  we  have  chosen  a  simple  two- 
dimensional  example  which  Is  based  on  a  passive 
ranging  experiment'*^. 

x(n)“jot(n-l)+u(n-l) , 

l(n)  ■h^  ( x(n)  ]+v(n) ,  (4) 


h  Ix(n)]'*tan“ 


j(n)-slnBn' 

(n)-cos6n 


where  the  transition  matrix  F  may  be  taken  as  the 
Identity  In  order  to  simulate  random  walk,  u  Is 
distributed  zero-mean  Gaussian  with  covariance  Q, 

V  Is  zero-mean  Gaussian  with  covariance  R  and  3  is 
the  angular  rotation  rate  of  the  sensor.  The 
mathematical  :!  (4)  arises  In  connection  with 
tracking  geometries  of  Figure  1  and  Figure  2.  In 
the  cose  of  Figure  1,  an  aircraft  sensor  S  Is 
orbiting  about  o  radarship  target  T  which  is  under¬ 
going  random  walk  motions  from  the  ocean  waves. 

In  this  case  the  initial  distribution  on  jt,  the 
position  of  the  ship,  is  taken  to  be  zero-mean 
Gaussian  with  identity  covariance  The  geometry 
of  Figure  1  more  closely  describes  the  situation 
studied  previously  than  that  of  Figure  2,  which 
Illustrates  the  more  important  (and  obviously  more 
difficult)  problem  of  locating  a  target  T  which  is 
separated  by  a  significant  distance  from  the  sen- 
cor's  circular  orbit.  Ihe  problem  described  by  (4) 


Is  In  either  case  one  of  locating  a  target  by  using 
only  bearing  sensor  Information  (passive  reception) . 
It  turns  out  that  the  previously  studied  problem  of 
Figure  1  Is  considerably  easier  than  originally 
suspected,  iind  a  new  experiment  comparing  all  pre¬ 
viously  proposed  methods  on  the  same  sample  tra¬ 
jectories  reveals  this  fact,  as  described  In  the 
next  section, 

5.  EXPERIME.NTAL  RESULTS 


In  order  to  adequately  study  the  passive  receiver 
problem  a  number  of  increasingly  difficult  cases 
must  be  considered.  The  first  case,  referred  to  as 
the  "old  problem,"  is  the  situation  illustrated  in 
Figure  1,  with  5=1.  Q*[o.05  O.ioj’ 

R«0,1,  and  N(0,I) .  The  second  problem, 
called  the  "new  problem,"  Is  identical  to  the  old 
problem,  with  the  important  exceptions  that 
^-N0],l),  R-0.01,  and  F-I.  Finally,  an 

example  with  a  multimodal  (l.e.  r.cngsusslan)  Irililal 
density,  described  In  Figures  3  and  4,  is  presented 
to  illustrate  what  happens  if  the  Gaussian  approxi¬ 
mation  is  radically  wrong,  r- - ; - r-r - 

Reproduced  from 
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It  turns  out  that  the  sensor  Is  allowed  to  have 


considerable  new  information  with  each  sample  if  it 
orbits  about  the  target  at  the  high  rate  of  S»1 
radian  per  sample.  Thus  it  is  not  surprising  to 
find  that  it  is  relatively  easy  to  design  an 
approximate  estimator  which  performs  very  close  to 
optimum.  A  straightforward  linearized  predictor 
performs  miserably,  it  happens,  because  of  the 
large  fluctuations  in  the  innovations  process 
Z|-h(x(n)  ]  as  the  sensor  completes  each  orbit.  In 
fact,  as  previously  observed'.^,  the  linearized 
filter  is  asymptotically  unstable  -  the  moan  square 
error  oscillates  with  ever  increasing  amplitude. 
This  fact  is  illustrated  in  Table  1,  where  in  the 
second  column  the  Monte  Carlo  average  of  the  sum 
squared  errors  for  100  sample  paths  is  presented 
for  an  iterated  linearized  predictor  wnich  is 
linearized  as  usual  about  the  current  filtered 


estimate  and  the  filter  update  is  iterated  until 
Che  sensor  linearizations  have  etabllized.  If  the 
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calculations  arc  continued  for  more  samples  in  each 
trajectory  then  the  instability  becomes  easily 
evident.  The  behavior  of  the  linearlaed  predictor 
is  so  poor,  in  fact,  that  it  is  vorsc  than  the 
"blind  predictor"  which  is  initialized  at  the  mean 
of  the  state  vector  and  ignores  the  data  entirely 
thereafter.  The  theoretical  mean-square  error 
performance  of  the  mean-state  predictor  Is 
illustrated  in  column  one  of  Table  1. 

If  tl^e  geometry  of  the  old  problem  is  exploited, 
it  is  possible  to  "fix  up"  the  linearized  pre¬ 
dictor  so  that  its  performance  la  almost  optimal. 
The  fixed  linearized  predictor  is  identical  to  the 
standard  linearized  predictor  with  the  addition  of 
logic  to  calculate  the  innovations  ^-h^fiCn)  ] 
modulo  2ii  in  the  range  The  sensational 

improvement  in  performance  of  the  fixed  predictor 
is  illustrated  in  column  three  of  Table  1. 

Finally  the  numerical  techniques  of  gaussian  sums®, 
and  floating  grids',  aie  included  iu  tlie  lust  two 
columns  of  Table  1.  If  allowances  for  the  fact 
that  the  one  standard  deviation  confidence  level 
in  a  Monte  Carlo  estimate  of  N  tracks  is  approxi¬ 
mately  20** /N,  where  Is  the  true  variance  and 
the  densities  arc  assumed  Gaussian,  then  for  100 
Monte  Carlos  the  resulting  variances  are  only  de¬ 
termined  within  one  standard  deviation  to  ±17.32 
percent  of  the  true  o^.  Thus,  to  within  the  exper- 
iointal  accuracy  available  in  100  Monte  Carlos, 
the  performances  illustrated  In  the  last  three 
columns  of  Table  1  are  essentially  equivalent. 

We  conclude,  therefore,  that  the  problem  is  not 
sufficiently  difficult  to  resolve  the  differences 
between  the  various  estimation  techniques. 

5.2  THE  NEW  PROBLEM 

r 

If  the  geometry  of  Figure  1  is  replaced  with  that 
of  Figure  2  then  it  no  longer  is  possible  to  fix 
up  the  performance  of  the  linearized  predictor  by 
a  simple  application  of  modular  arithmetic. 
Furthermore,  the  performance  of  the  numerical 
techniques  is  decidedly  better  than  the  linearized 
predictor  af  1 1  lustr.Ttcd  In  Table  2.  Tlie  Cau-islan 
sum  predictor  exi'crlenced  nr.  unexplained  transient 
error  in  tiu'  first  cstim.iti'  on  about  the  fiftieth 


Monte  Carlo  which  was  so  large  that  the  resulting 
Monte  Carlo  estimate  is  clearly  in  the  transient 
phase  in  the  Cable  for  the  first  estimate.  The 
recovery  is  stable,  however,  so  that  the  remaining 
values  in  column  three  are  reliable. 

5.3  MULTIMODAL  INITIAL  CONDITIONS 

In  the  previous  two  cases  the  relative  success  of 
the  llnearlzed-typc  predictor  can  be  related  to  the 
validity  of  the  Gaussian  approximation  to  the  con¬ 
ditional  density  functions.  The  simplest  way  to 
destroy  the  validity  of  the  Gaussian  assumption  is 
to  provide  a  nongausslan  initial  density  function. 
Consider,  for  example  the  detector  geometry  illus¬ 
trated  in  Figure  3.  If  there  were  a  sequence  of 
known  reflecting  ionospheric  layers  above ‘the 
aircraft  observer  and  we  were  given  an  apriori  dis¬ 
tribution  on  the  transmitter's  power,  then  it  is 
conceivable  that  we  night  want  to  integrate  over 
all  elevations  to  maximize  detectability,  thas 
Introducing  a  nulclmodal  range  ambiguity  as  shown 
in  the  figure.  Probabilistically,  the  initial 
condition  would  bo  obtained  by  taking  the  product 
of  the  multimodal  range  ambiguity  density  with  the 
bearing  ambiguity  density.  The  result  might  look 
as  in  Figure  4,  where  no  particular  scale  is 
intended.  As  soon  as  the  aircraft's  detector  cir¬ 
cuits  obtain  a  reliable  detection,  the  aircraft 
banks  left  into  a  circle  of  unit  radius  and  acti¬ 
vates  a  high  sensitivity  receiver  which  is  tuned  to 
one  elevation.  The  resulting  problem  parameters 

“rn  nen  'n  ““  ‘■[o.OOO  l.OOlj' 

Ho  ofs  0:05^' 

with  one  cross  range  standard  dcviatlc'n  correspond¬ 
ing  to  0.25  radian  in  beating  and  a  down  range 
standard  deviation  in  each  code  of  0.25.  If  the 
densities  arc  plotted  Isomctrically  as  viewed  from 
the  direction  shown  in  Figure  4,  then  some  of  the 
more  interesting  examples  from  the  movie  are  shown 
in  Figure  5,  with  sample  numbers  as  given.  As  seen 
in  Figure  5  the  density  first  peaks  up  on  the 
wrong  range  when  very  little  information  is  avall- 
.•<bl»‘  from  the  re.asurements ,  .ind  l.tter  makes  a 
di.in-atlc  chai  .;e  to  the  correct  mode  and  stabilizes. 
This  behavior  is  clearly  illustrated  in  the  solid 
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plot  of  F<«ure  6  of  the  resulting  prediction  error 
maghltude.  In  contrast  the  dotted  curve  In  Figure  6 
8h<MS  the  performance  for  the  same  sample  path  of 
the  linearized .predictor,  which  Is  almost  Identical 
to  the  mean  state  estimate  for  the  problem. 

Figure  6  Illustrates  two  Interesting  phenomena,  Ihe 
linearized  approximation  can  be  cade  useless  by  only 
a  large  uncertainty  In  initial  conditions.  On  the 
other  hand,  the  nonlinear  estimator,  which  contains 
the, more  detailed  Information  about  the  condi¬ 
tional  densities  may  still  perform  satisfactorily. 

In  addition,  the  optimum  systemcreates  a  certain 
amount  of  optimism  for  it.sclf  Initially  and  cal¬ 
culates  a  relatively  high  confidence  Incorrect 
estimate.  Experience  with  several  such  sample 
paths  has  resulted  In  the  conclusion  that  the 
situation  In  Figure  6  Is  just  about  the  worst  pos¬ 
sible  ease,  that  on  the  average  the  performance  la 
without  much  Improvement  for  about  half  the  itera¬ 
tions  <0.5  radian  of  sensor  lotatlons),  and  then 
proceeds  to  converge  at  about  the  some  rate  as  the 
previous  examples.  The  linearized  predictor,  on 
the  other  hand,  without  a  suitable  initial  condi¬ 
tion  doesn't  ever  converge  on  the  average. 

6.  CONCLUSIONS 

The  examples  have  shown  the  wealth  of  experimental 
experience  available  to  the  estimation  engineer  as 
a  result  of  numerical  studies  of  the  problem.  It 
is  possible,  for  example,  to  assess  the  validity  of 
given  analytical  approximations  and  the  value  of 
unused  Information  regarding  the  structure  of  the 
estimation  problem.  In  addition  It  Is  possible  to 
compute  practical  lower  bounds  on  performance  for 
nonlinear  eatlmators  In  some  small  dimensional 
problems,  a  fact  which  may  have  some  Important 
engineering  significance  to  designers  of  sensors, 
for  exainple,  who  must  know  the  theoretical  limita¬ 
tions  on  the  performance  of  a  given  design.  It  is 
anticipated  that  future  experiments  along  these 
lines  will  substantiate  this  conjecture.  In  the 
meantime  there  appears  to  be  many  possibr'  appli¬ 
cation  areas  to  explore. 
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Sample 

Kean-State 

Iterated 

Linearized 

"Fixed" 

Linearized 

Gaussian 

Sun 

Predictor 

Predictor 

Predictor 

Predictor 

1 

1.250 

4.334 

1.468 

0.757 

2 

1.388 

1.678 

0.974 

0.596 

3 

1.447 

•  0.651 

0.562 

0.457 

4 

1.537 

0.535 

0.550 

0.473 

5 

1.634 

0.828 

0.617 

0.518 

6 

1.734 

3.187 

0.619 

0.531 

7 

1.833 

4.927 

0.540 

0.471 

8 

1.933 

4.912 

0.618 

0.553 

9 

2.033 

4.248 

0.564 

0.578 

10 

2.133 

3.763 

0.619 

0.641 

Unstable 

Unstable 

Stable 

Stable 

Table  1.  Honte  Carlo  Averaged  Sum  Squared  Error 
Performance  for  Predictors  -  Old  Problem 


Floating 

Grid 

Predictor 

0.866 

0,674 

0.631 

0,408 

0.403 

0.464 

0.674 

0.454 

0.377 

0.443 

Stable 


Sas^le 

Mean-State 

Predictor 

Linearized 

Predictor 

Gaussian 

Sum  Predictor 

Floating  Grid 
Predictor 

1 

2.200 

1.665 

11.436* 

0.955 

2 

2.400 

1.928 

1.723 

1.020 

3 

2.600 

2.187 

1.495 

2.006 

4 

2.800 

2.229 

1.321 

1.083 

5 

3.000 

2.190 

1.208 

1.118 

6 

3.200 

2.145 

1.456 

1.359 

7 

3.400 

1.830 

1.625 

1.521 

8 

3.600 

2.105 

1.437 

1.370 

9 

3.800 

2.136 

1.423 

1.132 

10 

4.000 

2.349 

1.286 

1.316 

11 

4.200 

2.370 

1.431 

1.493 

12 

4.400 

2.114 

1.533 

1.515 

13 

4.600 

1.883 

1.475 

1.345 

14 

•4.800 

1.752 

1.398 

1.768 

15 

5,000 

1.752 

1.553 

1.441 

16 

5.200 

1.850 

1.617 

1.534 

17 

5.400 

1.627 

1.275 

1.222 

18 

5.600 

1.561 

1.271 

1.238 

19 

5.800 

1.620 

1.630 

1.314 

20 

6.000 

1.688 

1.803 

1.556 

*  transient 


Table  2,  Monte  Carlo  Averaged  Sun-Squared  Error 
Performance  for  Predictors  -  New  Problem 
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IONIZED  LAYERS 


Figure  3.  Multlmo'lal  Range  Anbigulty  - 
Detector  Gcocetry 


Figure  4.  Multimodal  Range  Arblgulcy  - 
Estimator  Ccococry 


Figure  6.  Absolute  E: 
.ind  I.lnearized  Prcdlc 


